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The  choice  of  an  optimal  design  for  parameter  estimation  in  a 

nonlinear  model  depends  on    the  values  of  the  model's  parameters.   This 

is  an  undesirable  feature  since  such  a  design  is  used  to  estimate  the 

same  parameter  values  which  are  needed  for  its  derivation.  To  overcome 

this  difficulty  a  method  is  developed  for  the  construction  of  a 

nonlinear  design  which  does  not  depend  on  the  nonlinear  model's 

parameters.   Such  a  design  is  called  a  parameter-free  design.   This 

method  is  an   extension  of  A.  I.  Khuri's  modified  D-optimality  criterion 

(Parameter-Free  Designs  for  Nonlinear  Models.  Technical  Report  No.  188, 

Department  of  Statistics,  University  of  Florida,  Gainesville,  Florida, 

1982).   It  is  based  on  using  a  proper  approximation  of  the  true  mean 

response  function  described  in  the  nonlinear  model  with  either  Lagrange 

or  spline  interpolating  polynomials.   Each  approximation  will  be  used  to 


derive  a  parameter-free  design.  Optimal  designs  obtained  on  the  basis 
of  these  interpolating  polynomials  depend  on  the  size  of  the  error 
associated  with  the  approximation  of  the  true  mean  response. 

A  method  for  the  construction  of  a  confidence  region  on  the  vector 
of  parameters  in  a  nonlinear  model  is  also  developed.  Unlike  most  other 
methods  which  are  available  for  that  purpose,  this  method  does  not 
require  specifying  initial  estimates  of  the  parameters.  The 
aforementioned  confidence  region  can  be  used  to  obtain  simultaneous 
confidence  intervals  on  the  nonlinear  model's  parameters.  These 
intervals,  however,  are  conservative  in  the  sense  that  their  joint 
confidence  coefficient  cannot  be  less  than  a  preset  value. 


CHAPTER  ONE 
INTRODUCTION 


1.1  Nonlinear  Models 


To  understand  what  we  mean  by  the  term  "nonlinear  model"  we  shall 
first  make  the  following  definition:  A  linear  model  is  the  one  in  which 
the  parameters,  the  quantities  to  be  estimated,  appear  linearly.  Such  a 
model  can  be  represented  by 

y  -   S0  +  81Z1  +  B2Z2  +  •••  +  3pZp  +  z       '         (I'D 

where  Sq,  S-^,  ...,  &Q  are  unknown  parameters,  the  Z;  are  any  functions 
of  the  input  or  independent  variables  x^,  ...,  x^ ,  y  is  the  dependent  or 
response  variable  and  e  is  an  unobservable  random  error. 

Any  model  which  is  not  linear  in  the  parameters  will  be  called  a 
nonlinear  model.  Examples  of  two  such  models  are 


y  =  exp[91  +  y2  +  e]  (1.2; 


y  =  q-4q-  [exp(-62x)  -  exp(-e1x)]  +  £   .        (1.3) 


The  models  (1.2)  and  (1.3)  are  both  nonlinear  in  the  parameters  81  and 
9?.  But  the  model  (1.2)  can  be  transformed  so  that  the  parameters 
appear  in  linear  fashion.  If  we  transform  the  model  by  taking  the 
logarithm  to  the  base  e  of  the  quantities  on  both  sides  of  the  equality 
sign,  we  obtain 


In  y  =  9.  +  9  x  +  e 


1.4) 


which  is  the  form  of  (1.1)  and  is  linear  in  the  parameters.  We  say  that 
model  (1.2)  is  an  intrinsically  linear  model  because  it  can  be  expressed 
in  the  linear  form  (1.1)  by  a  suitable  transformation. 

However,  it  is  impossible  to  transform  the  model  (1.3)  into  the 
linear  form  (1.1).  Such  a  model  is  said  to  be  intrinsically 
nonlinear.  Throughout  this  dissertation  all  models  mentioned  will  be 
intrinsically  nonlinear.  This  research  is  concerned  with  the  design  of 
experiments  for  parameter  estimation  and  the  construction  of  confidence 
regions  for  the  parameters  in  nonlinear  models. 


1 . 2  Design  of  Experiments  for  Nonlinear  Models 
Consider  the  nonlinear  model 


y  =  f(x;  9)  +  e 


1.5) 


where  y  is  the  measured  response,  f(x_;  _6_)  is  the  true  mean  response  at 
the  point  _x_,  where  _x_  is  the  vector  of  k  input  variables  (x-p  ...,  x^)1, 


_9_  is  the  vector  of  p  parameters  (9i,  8?,  ,..,    9Q)',  and  e  is  a  random 
error. 

When  there  are  n  observations  of  the  for;n  y^,  y?,  . . , ,  y_  with  y^ 
being  obtained  at  the  design  point  Xj ,  i  =  1,  Z,  . ..,  n,  then  we  can 
write  model  (1.5)  in  the  form 


y.  =  f  (x_.  ;  _8)  +  e. 


(1.6) 


where  Xj  =  (x-ji,  x-j?,  ...,  x^)'.  The  program  of  n  experiments  may  be 
represented  by  the  nxk  design  matrix  D_  =  [x^  •]•  The  j-th  member  of  the 
i-th  row  of  this  matrix  gives  the  level  of  the  j-th  variable  for  the 
i-th  experiment.  We  assume  that  the  random  errors  £i,  e?,  ...,  en  are 
independently  distributed  as  N(0,  ac) . 

The  design  criterion  for  parameter  estimation  in  a  nonlinear  model 
suggested  by  Chernoff  (1953)  is  to  choose  the  experimental  design  that 
minimizes  the  trace  of  [h_(x_>  0)  h_' (x_,  _§_)],  where  h_(x_,  _8)h_'  (x_,  Q)/a^  is 
Fisher's  information  matrix  of  order  pxp  and  h_'  (x_,  _9)  is  an  nxp 
derivative  matrix  with  ij-th  entry  3f (xj  ;  _e_)/88-.  The  problem  is  that 
the  Fisher  information  matrix  depends  on  the  unknown  parameters  _9_,  and 
one  cannot  determine  the  optimal  design  setting  for  estimating  _9  without 
knowing  the  actual  value  of  _6.  Chernoff  suggested  to  consider  locally 
optimal  designs  and  to  minimize  the  trace  of  [h_(x_,  _9q)  _h_'  (_x_,  ^q)Y  , 
where  _9q  is  the  value  of  _6  in  the  neighborhood  of  which  f  (xj  ;  _9_)  is 
assumed  to  be  approximately  linear  in  _9  for  i  =  1,  2,  ,..,  n.  Hamilton 
(1982)  pointed  out  that  this  criterion  is  not  invariant  under  changes  of 


scale  of  the  parameters.  Box  and  Lucas  (1959)  suggested  choosing  the 
experimental  design  to  maximize  the  determinant  |h_(x_,  _9q)  _hj  (x_>  On)  | 
with  respect  to  the  experimental  variable  settings  where  n,  the  number 
of  experiments,  is  equal  to  p,  the  number  of  parameters.  This  parameter 
estimation  criterion  was  called  D-optimality  by  Kiefer  and  Wolfowitz 
(1959).  Unlike  the  minimum  trace  criterion,  the  0-optimal  design  is 
invariant  to  scale  changes  of  the  parameters  (see  Cochran  (1973)  and 
Hamilton  (1982)). 

Atkinson  and  Hunter  (1968)  extended  the  Box  and  Lucas  results  to 
cases  where  n  is  greater  than  p.  In  some  cases  when  n  is  a  multiple  of 
p,  the  optimal  design  simply  results  in  n/p  replications  of  the  points 
which  would  be  given  by  the  n=p  Box  and  Lucas  approach.  However, 
Atkinson  and  Hunter  gave  a  counter  example  to  show  that  this  is  not 
always  true.  M.  J.  Box  (1968,  1970,  1971)  gave  a  number  of  additional 
results  particularly  on  the  subjects  of  replicated  design  points. 

Several  applications  and  theoretical  extensions  to  D-optimal 
designs  for  nonlinear  models  have  been  proposed.  Reviews  of  these  are 
given  by  Cochran  (1973)  and  St.  John  and  Draper  (1975). 

In  some  cases  it  may  be  true  that  initially  9i,  say,  will  be  known 
more  precisely  than  6?  at  the  outset  of  experimentation.  For  example, 
suppose  we  know  the  standard  deviations  of  their  prior  distributions  are 
one  and  ten,  respectively.  Draper  and  Hunter  (1967a  and  1967b) 
developed  a  procedure  for  the  design  of  experiments  for  parameter 
estimation  when  such  prior  information  is  available. 


For  the  nonlinear  model,  it  is  necessary  to  have  initial  estimates 
of  the  parameters  in  order  to  apply  the  previously  mentioned  design 
criterion.  If  these  initial  estimates  are  poor,  the  efficiency  of  the 
design  may  be  greatly  reduced  since  the  design  criterion  is  a  function 
of  the  parameter  values  that  are  employed.  A  more  efficient  approach 
would  be  to  treat  _6q  as  an  initial  estimate  in  a  sequential  procedure 
which  starts  by  selecting  a  p-point  optimal  design  on  the  basis  of  this 
estimate.  A  series  of  p  experiments  are  then  performed  at  the  optimal 
design  settings  and  an  estimate  of  _6  is  obtained.  Thereafter,  the 
experiments  can  be  designed  one  at  a  time,  using  the  current  best 
estimates  of  the  parameters,  which  can  be  recalculated  after  each 
experiment.  This  sequential  procedure  has  been  discussed  by  Box  and 
Hunter  (1965)  and  they  showed  how  to  add  points  one  at  a  time. 

M.  J.  Box  (1971)  added  sequentially  sets  of  p  points  each.  At  each 
stage  a  new  estimate  of  _9  is  obtained  and  used  to  generate  the  next  set 
of  design  points.  This  sequential  procedure  continues  until  the  values 
of  the  parameter  estimates  stabilize  along  with  the  values  of  the 
coordinates  of  the  design  points. 

The  sequential  procedure  described  above  can  help  improve  the 
efficiency  of  the  design  but  it  does  not  eliminate  the  dependency  of  the 
design  on  the  parameter's  values.  Furthermore,  Atkinson  and  Hunter 
(1968)  pointed  out  that  there  are  situations  in  which  it  is  not  feasible 
to  proceed  sequentially. 

Zacks  (1977)  considered  the  design  problem  in  a  3ayesian  framework 
which  is  constructed  so  that  the  expected  value  of  the  determinant  of 


Fisher's  information  matrix  with  respect  to  some  prior  distribution  of  _6 
is  maximized.  Bayes  sequential  designs  can  also  be  obtained.  In  this 
case  an  optimal  design  for  a  new  stage  of  experimentation  is  determined 
by  maximizing  the  expected  value  of  Fisher's  information  determinant 
with  respect  to  the  posterior  distribution  of  _9  given  the  results  from 
the  previous  stages.  The  Bayesian  approach  (sequential  and 
nonsequential)  eliminates  the  problem  of  dependency  of  the  optimal 
design  on  _9_,  but  it  introduces  dependency  on  the  prior  distribution. 

Hamilton  (1982)  introduced  a  design  criterion  based  on  the  second 
order  approximation  to  the  volume  of  the  confidence  region  of  the 
parameters  in  nonlinear  models.  This  criterion  is  more  efficient  than 
the  O-optimality  criterion  since  the  latter  is  based  on  the  first  order 
approximation  of  the  above-mentioned  volume.  However,  Hamilton's 
criterion  requires  good  initial  estimates  of  the  parameters  as  well  as 
of  the  error  variance. 

Khun'  (1982)  considered  two  types  of  designs  for  parameter 
estimation  which  do  not  depend  on  the  unknown  parameters  of  the 
nonlinear  model.  The  first  type  is  called  a  minimax  D-optimal  design. 
This  design  is  conservative  and  may  be  inefficient  as  pointed  out  by 
Khuri  .  The  second  type  of  design  is  more  efficient  than  the  minimax  D- 
optimal  one.  To  obtain  this  second  type,  an  adequate  approximation  to 
the  true  mean  response  with  a  Lagrange  interpolating  polynomial  was 
constructed.  This  approximation  was  used  to  derive  a  parameter-free 
design  on  the  basis  of  a  certain  modification  of  the  D-optimality 
criterion.   Khuri  noted  that  this  design  depends  on  the  size  of  the 


error  associated  with   the  approximation  of   the  true  mean   response.     Only 
a   nonlinear  model    with  a   single   input   variable  was   considered. 

In  this  dissertation  we  shall  extend  Khuri's  second  design  to  the 
cases  where  we  approximate  the  mean  response  function  with  a  Lagrange 
interpolating  polynomial  (for  the  case  of  two  input  variables)  and  with 
spline  interpolating  polynomials  (for  the  case  of  a  single  input 
variable) . 


1.3     Confidence  Regions   for  Parameters 
in   Nonlinear  Models 


Beale     (1960)     introduced     approximate     confidence     regions     for  the 

nonlinear    parameters    based    on    the    likelihood    ratio    criterion    when  the 

sample     size     is     large.        Williams      (1962)     presented     a     procedure  for 

obtaining   what    are    described    as    fiducial    intervals    and    regions    for  the 
nonlinear  model's   parameters.      He  considered  the  nonlinear  model 

y.    =   a   +   6f(xi ;    Y)    +   e.       , 


where  y  is  the  only  nonlinear  parameter.  Williams'  procedure  is 
applicable  for  a  model  with  only  one  parameter  appearing  nonlinearly. 
Halperin  (1963)  extended  Williams'  results  to  a  larger  class  of 
nonlinear  models  with  more  than  one  parameter  appearing  nonlinearly.  He 
showed  that  the  fiducial  intervals  and  regions  discussed  by  Williams  are 
also  confidence  intervals  and  regions.  Hartley  (1964)  suggested  a 
procedure  for  obtaining  an  exact  confidence  region  about  the  parameters 
in   a   general   nonlinear   model.     A  method   for   constructing   an 


approximation  of  the  true  mean  response  was  described.  This 
approximation  was  used  to  obtain  exact  confidence  regions  for  the 
parameters.  Sundararaj  (1978)  considered  a  method  for  the  construction 
of  joint  confidence  regions  on  a  vector  of  parameters  in  a  nonlinear 
model.  A  set  of  p  random  variables  was  constructed,  each  of  which  is  a 
least  squares  equation  for  estimating  the  parameters,  where  p  is  the 
number  of  parameters.  These  p  random  variables  were  used  to  construct 
the  joint  confidence  region  on  the  parameters. 

Several  methods  for  determining  confidence  regions  for  the 
parameters  of  nonlinear  models  have  been  proposed.  Reviews  of  these  are 
given  by  Gallant  (1976)  and  Shaw  and  Griffiths  (1979).  We  shall 
describe  briefly  some  of  the  methods  for  obtaining  the  confidence 
regions  on  the  parameters  in  nonlinear  models.  First  we  define  the 
error  sum  of  squares  for  the  nonlinear  model  given  in  (1.6)  as 


S(6)  =  z  (y,  -  f(xi:  £))' 

i=l 


(1.7) 


We  shall  denote  the  least  squares  estimate  of  8  as  9. 


1.3.1  The  Likelihood  Ratio  Method 


This  method  defines  the  confidence  region  for  9  by  the  expression 


S(9)-S( 


so: 


<  -2—  F(p ,  n-p ,  a) 

n-p   r 


which  provides  approximate  lOO(l-a)  percent  confidence  regions  for  8_, 
where  F(p,  n-p,  a)  is  obtained  from  a  table  of  the  F-di  stri  buti  on 
corresponding  to  a  significance  level  a  and  numerator  and  denominator 
degrees  of  freedom  p  and  n-p,  respectively  (see  Draper  and  Smith  (1931, 
p.  504)  and  Ratkowsky  (1983,  p.  30)). 

1.3.2  Asymptotic  Normality  of  the  Least  Squares  Estimator 

The  100(l-a)  percent  confidence  region  for  _6_  can  be  defined  as  (see 
Gallant  (1975a)) 


-6)1  C"1 


s(e: 


<  -2—  F(p,  n-p ,  a) 


n-p 


where      C    =    [h_(x_,     0_)   _h_'  (x_,    Q)Yl    and    h_(x_,  _9_)   h_' (x_,    Q)/a2    is    Fisher's 
information  matrix   as   defined   in   Section   1.2. 


1.3.3     Hartley's   Method 

Hartley  (1964)  developed  a  method  of  constructing  the  so-called 
"essentially  linear"  approximation  to  the  nonlinear  model  f {x_j  ;  _8)  of 
the  form 


P 
f(x    ;    e)    =      Z     w   (9)    U 

1  t=l    z 


where  wt(_9_)  are  continuous  functions  of  _9_  and  the  nxp  matrix  U  of  the 
u.jt  has  rank  p  and  does  not  depend  on  _9.  An  exact  100(1  -  a)  percent 
confidence  region  for  6  is  defined  as 


and 


10 


reg  SS/res  SS  <  -1—  F(p,  n-p,  a) 


reg  SS  =  (U'e)'  (U'U)"1(U'e 


res  SS  =  e'e  -  reg  SS 


1.4  The  Purpose  of  This  Research 

The  choice  of  an  optimal  design  for  parameter  estimation  in  a 
nonlinear  model  depends  on  the  values  of  the  model's  parameters.  The 
statistical  literature  suggests  several  procedures  for  obtaining  the 
designs  for  parameter  estimation  as  was  described  in  Section  1.2.  In 
most  cases,  it  is  necessary  to  have  initial  estimates  of  the  parameters 
in  order  to  apply  the  design  criterion.  If  the  estimates  upon  which  the 
design  is  based  are  poor,  the  design  may  be  inefficient. 

To  overcome  this  difficulty  we  present  a  design  criterion  for 
parameter  estimation  which  does  not  depend  on  the  nonlinear  model's 
parameters.  Our  design  criterion  for  parameter  estimation,  which  is  an 
extension  of  Khuri's  (1982)  procedure,  is  presented  in  Chapter  Two. 
Such  a  design  is  called  a  parameter-free  design.   The  approximations  to 


11 


the  true  mean  response  function  given  in  (1.5)  will  involve  a  Lagrange 
interpolating  polynomial  (for  the  case  of  two  input  variables)  and  the 
spline  interpolating  polynomials  (for  the  case  of  a  single  input 
variable)  in  order  to  obtain  these  parameter-free  designs.  Chapter 
Three  contains  a  method  for  the  construction  of  confidence  regions  on 
the  parameters  in  nonlinear  models.  This  method  does  not  require 
specifying  initial  estimates  of  the  parameters.  Using  the  confidence 
intervals  on  the  mean  response  functions  in  conjunction  with 
Bonferroni's  inequality,  a  conservative  confidence  region  on  the 
nonlinear  model's  parameters  is  constructed.  In  the  final  chapter,  we 
make  some  concluding  comments  concerning  our  methods. 


CHAPTER  TWO 
DESIGNS  FOR  NONLINEAR  MODELS 


2.1  Introduction 


In  this  chapter  we  introduce  techniques  for  constructing  the 
polynomial  approximations  of  a  nonlinear  model.  Each  approximation  will 
be  used  to  obtain  a  design  for  parameter  estimation  which  does  not 
depend  on  the  parameter  vector  _6.  We  shall  restrict  our  consideration 
to  the  cases  of  one  and  two  input  variables. 

Consider  the  nonlinear  model 

y  =  f(x_;  _9)  +  e   ,  (2.1) 

where  x_  =  (x^,  x?,  ...,  xk)'  is  a  vector  of  k  input  variables,  _6_  = 
(81,  9jp,  ...,  9  )'  is  a  vector  of  p  unknown  parameters,  and  e  is  a 
random  error.  To  estimate  the  elements  of  8,  a  series  of  n  experiments 
(n>p)  is  performed.  In  each,  a  value  of  the  response  is  observed  at  a 
particular  setting  of  _x_.  The  program  of  n  such  settings  is  represented 
by  the  n  x  k  design  matrix,  D_  =  [x^],  i  =  1,  2,  ...,  n;  j  =  1,  2,  ..., 
k.  We  assume  that  the  random  errors  e-t,  e?,  ,..,  e  are  independently 
distributed  as  N(0,  a2) . 
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We  now  discuss  the  problem  of  finding  the  polynomial  approximations 
to  the  given  nonlinear  models  of  one  and  two  input  variables. 

2.2  Approximation  and  Interpolation  of  Functions 
One  of  the  simplest  ways  to  approximate  a  function  defined  on  an 
interval,  is  to  obtain  a  polynomial  which  takes  the  same  values  as  the 
function  at  a  certain  number  of  points  in  the  domain  of  the  function. 
This  procedure  is  called  interpolation.  Szidarovszky  and  Yakowitz 
(1978)  pointed  out  that  polynomials  provide  a  simple  way  to  approximate 
or  interpolate  functions.  By  way  of  illustration,  let  us  state  a  well- 
known  theorem  of  Weierstrass. 

Theorem  2.1  (The  Weierstrass  Approximation  Theorem) 

Given  any  interval  [a,b],  any  real  number  6  >  0,  and  any  real - 
valued  continuous  function  f(x)  on  [a,b],  there  exists  a  polynomial  p(x) 
such  that 

|f(x)  -  p  (x)  |  <  6   for  all  xe[a,b] 


This  result  may  be  proved  in  a  number  of  ways:  see,  for  example,  Davis 
(1963)  and  Rivlin  (1969). 

In  this  dissertation  we  shall  center  our  attention  on  the 
approximation  of  a  nonlinear  model  by  considering  a  Lagrange 
interpolating  polynomial  (for  the  case  of  two  input  variables)  and 
consider  spline  interpolating  polynomials  (for  the  case  of  a  single 
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input  variable).  Each  approximation  will  be  used  to  obtain  a  parameter- 
free  design  on  the  basis  of  a  modification  of  the  O-optimality  criterion 
which  was   proposed   by   Khuri    (1982). 

2.2.1     Lagrange   Interpolation 

One-Dimensional    Problem 

The  Lagrange  polynomials  are  among  the  simplest  and  most  practical 
of  the  interpolating  polynomials.  Given  m+1  distinct  real  numbers 
a<yn<y,<. .  .<y  <b  and  a  function  f(x)  defined  at  yQ,  y^,  ...,  y^,  let 
g(x)   be  a  polynomial    of  degree  m  or  less   for  which 


g(yi)  =  f(yi)   ,  1  =  o,  l,  2,  ...,  m 


(2.2) 


This  polynomial  can  be  written  in  Lagrange  form 


g(x)  =  I    f(yi)i,(x) 

i=o   n  n 


2.3] 


where 


-1 


l-(x)  =  n  (x  -  y.)(y.  -  y.)   ,  i  =  0,  1, 
j»0      J       J 
j*1 


,  m. 


(2.4) 


is  the  i -th  Lagrange  polynomial  of  degree  m  (see,  for  example,  Young  and 
Gregory  (1972,  Ch .  6),  and  deBoor  (1978,  Ch .  1  and  2)).  We  note  that 
Ij(yj)  =  0  if  i  *  j  and  1  i  (yi )  =  1.   The  points  yQ,  y1?  ...,  ym    are 
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called  interpolating  points  and  the  polynomial  g(x)  is  called  Lagrange 
interpolating  polynomial  of  degree  m.  Such  a  polynomial  always  exists 
and  is  unique  as  can  be  stated  in  the  following  theorem. 

Theorem  2.2 

There  exists  a  unique  polynomial  g(x)  of  the  form  (2.3)  for  which 

g(y0)  =  f(*o)'  9(yi)  -  f(yi).  •••,  g(ym)  =  f(ym). 

The  proof  of  this  theorem  can  be  found  in  Prenter  (1975),  and 
Szidarovszky  and  Yakowitz  (1978). 

The  error  resulting  from  the  use  of  Lagrange  interpolation  is,  by 
definition,  the  difference  between  f(x)  and  g(x),  where  g(x)  is  the 
Lagrange  polynomial  of  degree  m  interpolating  f(x)  at  m  +  1  distinct 
real  numbers  yn,  y^,  ...,  y_.  We  assume  that  f(x)  has  continuous 
derivatives  up  to  order  m  +  1  with  respect  to  x  over  the  interval 
a<x<b.  Then  the  interpolation  error  is  given  by 


e(x)  =  f(x! 


g(x; 


m 

n  (x 

1*0 


yi)f(rT1+1)(?)/(m+l)! 


:2.5) 


where  a  <  ?  <  b,  xe[a,b]  and  f^+1^(x)  =  3m+1f (x)/3xm+1  (see,  for 
example,  Prenter  (1975,  p.  36)  or  deBoor  (1978,  p.  8)).  As  5  is  not 
explicitly  known,  formula  (2.5)  cannot  be  used  for  the  calculation  of  an 
exact  value  of  e(x)  but,  as  we  shall  now  see,  it  can  be  used  for 
obtaining  an  upper  bound  for  the  absolute  value  of  e(x).  If  we  let 


w(y0.  yv 


,  y  )  =  max  I  n  (x  -  y . 
'    a<x<b  i=0 
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an  upper  bound   for    |e(x)|    is   given   by 

max    |e(x)|    <  w(yn,  y. ,    ...,  y   )   max    \f'n]+     (x)  |/(m+l; 
a<x<b  u       L  m  a<x<b 


'2.6} 


Notice  that  the  interpolation  points  influence  the  size  of  the  upper 
bound     on     the    error     through     the     function     w(yn,     y^,     ...,    y   ).        The 

function   w(yn,   y-i.    >   ym)    can   be   made   small    by    suitable  choice   of  the 

points  y.j .  In  fact  the  values  of  y^  which  minimize  the  function 
w(y{]>  yi .  •••»  ym)  can  be  obtained  by  using  the  following  result  given 
in  de3oor  (1978,  p.  26-31):  the  function  w(yQ,  y1$  ...,  y  )  in  (2.6) 
attains  a  minimum  value  of  2(b  -  a)m  /4m+1  if  the  interpolation  points 
are  chosen  as  the  zeros  of  the  Chebyshev  polynomial  of  the  first  kind 
and  of  degree  m+1  for  the  interval  a<x<b.  These  zeros  are  given  by  the 
formul a 


a+b       ra-th       r,-       1\   n  i     .   „     , 
yi    =  ~  "   C— 2— )cos  L(  "I    +   2-feJ'    1=0'    l'    •'•,   m        '  {-ZJ' 


From    (2.6)    and   the   fact    that   the   minimum   value   of  w(yg,  y-,,    ...,  y   )    is 
2(b   -  a)m+1/4m+1  we   obtain 


max    |e(x)|    <   2(b   -  af+l  max    |  f (m  +  1)  (x)  |  /[4m+1  (m  +  1) !  ]      .       (2.8) 
a<x<b  a<x<b 
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Two-Dimensional  Problem 

We  now  look  at  the  problem  of  approximating  functions  of  two 
variables  with  Lagrange  polynomials  in  two  variables.  Suppose  that  we 
are   given  a  rectangular  region 


Rx  =  {(xr  X?J:   al  *  xl  <  br  a2  <  x2  *  b2} 


a,  <  yn<. .  .<  y  <b.   and  n 
x,   1  JQ  Jm.   1       x? 


a2  <  zQ  <  Zj  < 


mp    ■:: 


Let   n 

x 

be  partitions  of  the  intervals  [a^,b^]  and  [ap, bo],  respectively,  where 
yQ»  y\,  •••>  ym  and  Zq,  z^,  ...,  zm  are  m^  +  1  and  m2+l  distinct  real 
numbers.  The  set  of  points 


n  =  {(yi  ,  z.)  :  0<i<mi,  0<j<m2} 


partitions  R  into  a  family  of  subrectangl es 


Rij  :  *1-1  <  xl  <  V  zj-l  <  x2  <  zj 


J  =  1.  ••  •,  m^ ;  j  =  1,  .. .,  m 


y0  =al'  \  =  bl'  Z0  =  V  Zm2  "  ^ 


Let  f(xj,  x2)  be  a  function  defined  on  R   and  let  g(xi,  x?)  be  a 
polynomial  of  degree  m^  in  X]_  and  m2  in  X2  such  that 
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g(yr  z.)  =  f(yr  z.)  .  0<i<mi;  0<j<mn   .       (2.9; 


This  polynomial  can  be  written  in  Lagrange  form 


g(x, ,  x?)  =  I        f(y.  ,  z.)  1,.-(x15  x  ' 

1   c  0<i<m1       J    J 

0<j<mo 
where 


•2.101 


1ij^xl'   x2^   =  ]i  ( x i ) "•  j  ( x 2 ' 


'2.111 


]i{xi]  \t0  (xi  -yk)(yi  -v1  ■  i=°> i 

k*i 


lj(x2)    ^(x2-zu)(2.    -zj"1,   J-0.    1. 


'1       ' 


2.121 


l-j(x^)  and  1  ,-  ( x  2 )  are  the  i,  j-th  Lagrange  polynomials  of  degree  m^  and 
mo,  respectively  (see  Prenter  (1975,  p.  119)).  We  note  here  that 


1ij^k'  zu> 


'  1     if  i  =k  and  j=u  , 
0     otherwise. 
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Theorem   2.3 

There  exists  a  unique  polynomial  of  degree  mi  or  less  in  xi  and  m? 
or  less   in   x2  which   satisfies    (2.9). 

The  proof  of  this  theorem  can  be  found  in  Prenter  (1975,  p.  120). 
We  now  look   at  the  error  of  interpolation.     Let  us   assume  that   f(x-.,   x2) 

has   continuous   derivatives   up   to   order  mj+1   and  m2+l   with    respect   to  x, 

2 
and    x2,    respectively,    over    the    rectangular    region      R    ;      then    an    upper 

bound   for  the  absolute  value  of  e(x^,   Xo)    is   given   by 


Wbi 

a2<x2<b2 


|e(x1}  x2, 


a, <x, <b . 
a2<x2<b„ 


f(xj,  x2)   -  g(xr  x2! 


\^o'  ■•■■  ^ 


max        |f  (x, ,   x?) 

a1<x1<b1     xl  L 

a?<x„<b2 


,   z 


(n>2+l)! 


(m2+l) 

max        |f  (x,  ,    x,) 

al<xl<bl 
a?<x,,<b,, 


T   A2* 


(v0, 


»ym  K 


'0' 


,z 


(mi+l)!(m2+l)! 


max 


(n^+1)  (m2+l) 

fy  (X-,  ,     X?)f  (x,  ,     X-)| 

a1<x1<b1  xl  i       c     x2  x       ^ 

a?<x2<b? 

(2.13) 
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where 


(m  +1) 
fl       (xr  x2) 


3  f(Xp    Xg) 


3x 


171,+! 


(m2+i; 

fx  (xis  x9) 


r  *2; 


m„+l 

3     f(xr  x2: 


3x 


m2+l 


and 


"x^O' 


.  y„ 


a ,  <x .  <b .   i=0 


-Q,      • 


,    2 


max 
a2<x2<b2  j=0 


? 
IT  (x, 


z  . 

J 


(see  Prenter  (1975,  p.  125)). 

Consider  now  the  upper  bound  of  the  error  in  (2.13).  Suppose 
f(xj,  x2)  is  a  function  of  xj  alone,  that  is  f(xi,  x?)  is  reduced  to  a 
one-dimensional  function.   Then  f  ^        (x,,  x5)  =  0,  0<j<m„.  Thus,  the 

X2       1    <J  2 

upper  bound  in  (2.13)  is  reduced  to 


(nii+1) 

max   |f      (x , 
a1<x1<b1   1 

(m^l)! 


"x^O- 


.  yn 


which  is  the  same  as  the  upper  bound  in  (2.6). 

Now  if  yQ,  y,,  ...,  ym  are  chosen  as  the  zeros  of  the  Chebyshev 
polynomial  of  the  first  kind  and  of  degree  iiii+l  for  the  interval 
[a^,bj],  and  Zq,  z^,  ...,  zm  are  the  zeros  of  the  Chebyshev  polynomial 
of  the  first  kind  and  of  degree  m2+l  for   the  interval  [a2,b?],  that  is, 
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>',     :  -^~  -   (—g— )cos[(i    +  ?)  m^'   1=0'    l«    •••-  mi     • 


i 


.^.(^costtj.l,.^.],  j-o.  1, 


(2.14) 


then    the    functions    wx   (yQ,     ...,    y     )    and    wx   (zQ,    ...,    z     )    in     (2.13) 


mi+1     mi+1 
1     //i   i 


Olp  +  1       ITlp+1 


attain  minimum  values  of  2 ( b -^  -  a^)    /4    and  2(b?  -  a?)  '  /4 
respectively.    In  this  case,  the  upper  bound  of  the  error  of 


interpolation  becomes 


(■n,+l) 
2(b1   -  aL)     l  (m1+l) 

max        \e(xv  x2)|    < — max        |f 

al<Xl<bl  (m1  +  l)!4   l  Wbl        l 


a2<x2<b2 


1 


a?<x?<b? 


, x  ,  ,    Xp , 


m?+l 
2(b2-a2)  (m2+l) 

m7+T        ™*        lfx9  (xl'   x2) 

(m2+l)!42    Vxl<bl   2 
a2<x2<b2 


m,  +1  r 

Mbj-aj)         (b2-a2) 


it,  +1  m?+l 

(mj+1)  (m2+l) 

m.^+T      fJV  (V  x2)fx,  (xl'  x2; 

(m1+l)!(m2+l)!4  !     2        al<xl<bl     X  2 

a?<x?<b2 


2.15) 
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We  have  described  the  Lagrange  interpolating  polynomial  for  the 
cases  of  one  and  two  input  variables.  The  construction  of  this 
approximation  is  fairly  simple.  However,  a  drawback  to  this 
approximation  is  that  it  is  not  smooth,  since  it  cannot  be 
differentiated  at  the  interpolating  points.  In  the  next  section  we 
introduce  a  smoother  approximation.  This  approximation  is  called  a 
spline  interpolating  polynomial. 


2.2.2  Spline  Interpolation 

In  the  last  30  years  or  so  following  the  pioneering  work  of 
Schoenberg  (1946),  piecewise  polynomial  functions  have  become  an 
important  theoretical  tool  for  the  approximation  of  functions,  since  in 
many  cases  the  behavior  of  a  function  in  one  region  may  be  totally 
different  from  its  behavior  in  another  region  and  one  polynomial  may  not 
be  appropriate  to  approximate  that  function.  The  most  popular  class  of 
approximating  functions  is  the  class  of  spline  functions.  Therefore, 
instead  of  trying  to  approximate  a  function  (of  possibly  several 
variables)  over  the  entire  region  by  one  polynomial  of  a  certain  degree, 
one  partitions  the  region  into  several  subregions,  and  approximates  the 
function  on  each  of  these  subregions  by  (possibly)  a  different 
polynomial.  One  of  the  important  ideas  in  using  this  process  is  the 
method  in  which  the  region  is  subdivided.  In  one  dimension,  the  method 
consists  of  defining  points  (called  knots)  whose  positions  generally  are 
not  known  beforehand.  Ideally,  we  would  like  to  have  these  knots  as 
variables  that  can  be  shifted  into  the  optimal  positions  as  the 
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computation  proceeds.  However,  this  turns  out  to  be  a  difficult 
nonlinear  problem  and  so  the  knot  positions  are  usually  fixed  in 
advance.  Throughout  this  dissertation  we  assume  that  the  knot  positions 
are  preselected,  so  that  the  problem  is  a  linear  approximation  problem; 
furthermore,  for  simplicity  we  shall  restrict  our  consideration  to  the 
case  of  equidistant  knots.  Further  discussion  of  the  positioning  of  the 
knots  can  be  found  in  Powell  (1967),  Rice  (1969)  and  Wold  (1974). 
Estimation  of  the  knot  locations  is  discussed  by  Gallant  and  Fuller 
(1973). 

We  shall  now  discuss  the  spline  interpolating  polynomial  of  degrees 
1,  2  and  3  (linear,  quadratic  and  cubic  splines).  We  shall  consider 
only  the  case  of  a  single  input  variable. 

We  consider  the  interval  [a,b]  and  divide  it  into  h+1  pieces  using 

h  points  or  "knots"  x.  ,  x„,  ...,  x,   such  that   a  =  xn<x1  <• .  ,<x.  < 
1   2       h  ij  l      h 

x,.  =  b.  Generally,  a  function  which  is  equal  to  a  polynomial  of 
degree  at  most  m  on  each  subinterval  (tj,t--+i),  i  =  0,  1,  ...,  h,  can  be 
written  as 


m    .    m  .         m 

g(x)  =  Z   a.x1  +  £  a  .  (x  -  t   )]    +  ...  +  I   a,  .  (x  -  x  )]    , 
i=0      i=0  n      l  i=0  ni      n 


2.16) 


where 


(x-x  )+  =  (x-i  .)'  ,    X  >  T  . 


1 

J'   '  "   J 


(i  =  1,  2, 
=  0        ,   x  <  x. 
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A  function  g(x)  is  called  a  spline  of  degree  m  if  it  is  a  polynomial  of 
degree  (at  most)  m  on  each  of  the  h+1  intervals  (t.,t.  .),  i  =  0,  1, 
...,  h,  and  has  m-1  continuous  derivatives  at  x.- ,  i=l,  2,  ...,  h.  By 
imposing  the  restriction  that  g(x)  in  (2.16)  has  m-1  continuous 
derivatives  at  i|,  i  =1,  2,  ...,  h,  we  obtain  a  spline  of  degree  m  of 
the  form 


m     .    h 
g(x)  =  I     a.x1  +  I     b. 
i=0  7     i=l  ] 


(2.17; 


(see,    for  example,    Studden  and  Van   Annan    (1969),    Studden    (1971)   and   Park 
(1978)). 

For    illustration,    suppose   cubic-cubic   polynomials    are   adequate    for 
a   response   function   f  given   in    (2.1)   over  the   interval      a<x<b. 


2  3 

g(x)    =  a10   +  a^x   +  a1?x      +  a^x        ,   a    <   x   <• 


2  3 

a__   +  a_.x    +   cu_x      +  a~,x         ,    x   <    x    <  b 


'2.13) 


We  can  write  g(x)  in  (2.18)  as  the  single  equation  model 


2      3         0 
g(x)  =  aQ  +  a^  +  a2x  +  a3x  +  a4(x-x)+  +  a5(x  -  t)  + 


+  ag(x  -  t)2+   +  ay(x  -  x)l       , 


'2.19! 
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where 


a0  =  a10'  al  =  all'  a2  =  a12'  a3  =  a13 


a4  =  a20  "  a10  +  x(a21  "  all>  +  x2(a22  "  a12)  +  x3(a23  "  a13; 


a5  =  a21  "  all  +  2x(a22  '  a12)  +  3x  (ct23  "  a13: 


a6  =  a22  "  a12  +  3x(a23  "  a13] 


We  assume  that  g(x)  in  (2.19)  is  continuous  and  has  first  and 
second  continuous  derivatives  at  x=t;  thus,  from  the  condition  of 
continuity  at  x=x,  we  can  obtain  from  (2.18)  that 

a10  +  allx  +  a12x2  +  a13x3  =  a20  +  a21x  +  a22x2  +  a23x3  ' 


or 


°20   "  a10  +  x(a21   "  all>    +  x2(a22   '  a12)    +  x3(a23   '  a13)    =  °        ' 


which    implies    a^  =   0  in    (2.19).      From  the   condition  of   first   continuous 
derivative   at   x=x,  we   have   from   (2.18)    that 
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a21  -  au  +  2x(a22  -  a^)  +  3x  (a^  -  o^)  =  0 


hence  a^  =  0  in  (2.19).   Finally,  we  impose  that  g(x)  in  (2.19)  has 
second  continuous  derivative  at  x=t,  and  we  obtain 


a22  "  a12  +  3x(a23  "  a13}  =  °   • 


which  implies  that  a,-  =  0.  Thus  g(x)  in  (2.19)  becomes 


2     3  3 

g(x)  =  a„  +  a,x  +  a~x  +  a^x  +  a7(x  -  x) 


Throughout  this  dissertation  all  the  splines  of  degree  m  have  m-1 
continuous  derivatives  at  the  knots.  A  spline  function  of  degree  m  with 
h  interior  knots  is  represented,  in  general,  by  a  (possibly)  different 
polynomial  in  each  of  the  h+1  intervals  into  which  the  h  interior  knots 
divide  the  real  line.  As  each  polynomial  involves  m+1  parameters,  the 
spline  function  involves  a  total  of  (m+1) (h+1)  parameters.  However,  the 
continuity  conditions  impose  certain  constraints  on  these  parameters. 
At  each  interior  knot,  t^  ,  i=l,  2,  ...,  h,  the  spline  function  has 
continuous  derivatives  of  orders  0,  1,  ...,  m-1.  Thus,  m  constraints 
are  imposed  at  each  knot.  As  there  are  h  interior  knots,  a  total  of  mh 
constraints  are  imposed.  The  number  of  free  parameters  involved  in  the 
spline  function  is  the  total  number  of  parameters  minus  the  number  of 
constraints,  that  is,  (m+l)(h+l)  -  mh  =  m+h+1.   We  note  here  that  the 
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number     of 

free     parameters     is     equal     to     th 

e     numbe 

r     of 

a  's 

and 

b's 

in    (2.17). 

2.2.2.1     L 

i  near  Spl i  ne 

Let 

a   =  V-r^t^. 

..<Vxh+1=b       be 

partition     pc 

ints 

for 

the 

interval     [ 

a,b]    and    let 

f(x)     be    a    function 

def i  ne 

i    at 

T0' 

Tl' 

. . ., 

Th+1-       We 

construct    a 

linear    spline    function    as 

follows . 

On 

?ach 

interval 

I*t.  ,t.      ],      we 
-  i      i+lJ' 

choose    the    interpo' 

ant    g 

:o    f 

to    consist 

of 

1 inear  pieces, 

g(x 

)    =   P.(x)    =  a. 

+    V>     T.      <X     <     T.+1, 

i=0,    1 

,    ... 

,   n 

•    (2 

.20) 

The     i -th 

linear     piece 

Pi      is     made     to     satisfy 

the 

interpol atory 

constraints 

w    = 

f(t.)             ,   i=0,    1, 

2,    ..., 

h, 

ph^h+l) = 

f<W          ' 

(2 

21) 

VW- 

Pi+1(-+1)    ,    1-0.   1, 

2,    ..., 

h-1 

(see     Ahlberg,     Ni 1  son     and     Walsh     (1967,     p. 

109)). 

The 

resul t 

i  ng 

piecewi  se 

inear   function 

g(x)    agrees   with   f(x 

)       at      Tg 

'    Tl' 

Th  +  1 

and 

is    continuous    at    t,,    t?, 

...,    Th.       By    using 

(2.17) 

with 

m  =  l, 

g(x) 

in 

(2.20)   can 

be  wri  tten    in 

the    form 
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g(x)  =  aQ  +  ajX  +  b^x  -  ^)  +  +  b2(x  -  x?)  +   +  ...  +  bh(x 


(2.22) 
where  9q,  a-,,  b-^,  ...,  b^  can  be  obtained  by  using  the  conditions  in 
(2.21).  We  note  that  the  number  of  parameters  in  (2.22)  is  equal  to  the 
number  of  free  parameters  which  is  h+2. 

Consider  now  the  error  resulting  from  the  use  of  a  linear  spline. 
We  assume  that  the  function  f(x)  has  continuous  derivatives  up  to  order 
2  with  respect  to  x  over  the  interval  a<x<b.  Then  an  upper  bound  for 
the  error  is  given  by 


max  |g(x; 
a<x<b 


f(x)|  < 


, At)  max 
a<x<b 


f(2)(x: 


/8 


2.23] 


where 


At  =  t. , .  -  t. 
i  +1    i 


i=0,  1 


and 


f(2)(x)  =  32f(x)/3x2 


(see  deBoor  (1978,  p.  40)). 


2.2.2.2  Quadratic  Spline 

Let  f(x)  be  a  function  defined  at  in,  t.  ,  ...,  t,  .   with  a  =  t„< 

U   1       h+1  0 

t  <...<x  <Th   =  b.   We  choose  the  interpolant  g  to  f  to  consist  of 
quadratic  pieces, 
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g(x)   =  Pi  (x)   =  o^    +  67-x  +  y^x     ,   x.    <  x  <  ii+1,   1=0,   1.    ...,  h     . 

(2.24) 

The     i -th     quadratic     piece     Pi      is     made     to     satisfy     the     interpolant 
constraints 


P,(t, 


^Ul)       =  WW    ■    1-0.   1.    ....   n-1 


where 


Pi  '<w  =  WW  >i=0>  l  ■•••  h-! 


P0(XQ)     -f(T0)         , 


vw  =  f<W 


P^^O    =  f(^-)        ,       1-0,   1,    ...,   h. 


Ci    ■   (^    +  Ti+1)/2,   1*0,    1,   2,    ...,   h,    and 


P       (Ti+i)   =  C3Pi(x)/3x]y_T        ,   1=0,    1,    ...,   h-1 
*     i+1 


2.25) 


(see  Ahlberg,  Nilson  and  Walsh  (1957,  p.  109)  and  deBoor  (1978)).  The 
piecewise    quadratic     function     g(x)     interpolates     the    function     f(x)     at 

tq»  Tn+i>  ?o'  ?1'  '"'  S  and  ^as  conf-inuous  derivatives  of  orders  0 
and  1  at  Xp  x2>  ....  x^.  By  using  (2.17)  with  m=2,  the  function  g(x) 
in    (2.24)   can   be  written   as 
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g(x! 


2  2  2 

aQ  +  a,x  +  a^x  +  b1(x-t1)+  +  b2(x-x2)+  + 


b  fx  -  x  ) 

V     h;  + 


(2.26) 
where  a0,  a,,  a2,  bj,  b?,  ...,  bh  can  be  obtained  by  using  the 
conditions  in  (2.25).  The  number  of  free  parameters  is  equal  to  the 
number  of  parameters  in  (2.26)  which  is  h+3. 

To  obtain  the  error  |g(x)  -  f(x)|,  we  assume  that  the  function  f(x) 
has  continuous  derivatives  up  to  order  3  with  respect  to  x  over  the 
interval  a<x<b;  then  the  upper  bound  of  the  error  resulting  from  the 
use  of  quadratic  spline  is  given  by 


where 


and 


max    |g(x)    -  f(x)|    <   (At)3  max    |f(3)(x)|/8 
a<x<b  a<x<b 


At  =  t.    .    -  t.        ,  i=0,   1,   2,    ...,  h 


f(3)(x)    =  33f(x)/3x3 


'2.27) 


(see  deBoor   (1978,   p.   81) 


2.2.2.3     Cubic   Spline 


Let    a   =  Tq   <   t,    < 


.     <      T,       <      T 


l   <    x^+^    =   b   be   partition   points    for   the 


interval     [a,b]    and     let     f(x)     be    a     function    defined    at     Tq,    t,,     ..., 
THi  +  i«     We   construct   a  piecewise  cubic   interpolant   g  to   f  as   follows.     On 
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each  interval  [t^,t^+,],  let  g  agree  with  some  cubic  polynomial  of  the 
form 


2      3 
g(x)  =  P.(x)  =  a.  +  0.x  +  y-x  +  6.x  ,  x.  <  x  <  x.+1,  i=0,  1,  ...,  h 


The   i-th   cubic  piece   P.    is  made  to   satisfy  the   conditions 


2.28) 


PAt 


,(t1+1)        =   WW'    i=0'    L    —    ^ 


P{1}(^+1)   =?\\\^^>   i=0>   L    —   ^ 


p{2)<W  -P{?}(t1+1),  1=0,  1,   ....  h-1 


P.(x.)    =  f(x.),    i=0,    1,    ....    h 


:2.29! 


Vvi'  =  f<W    • 


o    l  o;  l  c 


h         <W      =     f  "hrl 


where 


(J) 


(W    =    [3JPi(x)/3xj]  ,    j    =    1,    Z, 

11  ]  i+1 
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and 


f(1)(t.)    =    [3f(x)/3x]x    =  x 


(see  Ahlberg,  Ni 1  son  and  Walsh  (1967,  p.  109)  and  Prenter  (1975, 
p.  78)).  The  piecewise  cubic  function  g(x)  interpolates  the  function 
f(x)  at  Tq,  T-^,  ...,  xh  +  ^  and  has  continuous  derivatives  of  orders  0,  1, 
and  2  at  t^,  t^,  ...,  t^  .  By  using  (2.17)  with  m=3,  the  function  g(x) 
in    (2.28)   can  be  expressed  as 


g(x)   =  aQ  +  axx   +  a2x     +  a3x°  +  b^x   -  t^J  +  b2(x   -  tJ* 


'2.30) 


.+bh(x    -xh)J 


where  3q,  a^,  a?,  93,  b^,  b?,  ...,  bh  can  be  obtained  by  solving  the 
conditions  in  (2.29).  The  number  of  free  parameters  is  equal  to  the 
number  of  parameters  in  (2.30)  which  is  h+4. 

We  assume  that  f(x)  has  continuous  derivatives  up  to  order  4  with 
respect  to  x  over  the  interval  a<x<b.  An  upper  bound  on  the  error 
resulting  from  the  use  of  cubic  spline  is  given  by 


max  I g ( x ) 
a<x<b 


f(x)|  <  5(Ax) 


max  |f 
a<x<b 


!x)|/384 


(2.31; 


where 


Ax  =  T 


i+1 


,  i=0,  1,  2,  ...,  h, 
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and 


f^(x)  -  34f(x)/3x4 


(see  Hall  (1968)  and  deBoor  (1978,  p.  68)). 

2 . 3  Optimal  Designs 
Khun'  (1982)  obtained  an  adequate  approximation  to  the  mean 
response  function  with  a  Lagrange  interpolation  polynomial  and  used  this 
approximation  to  derive  a  parameter-free  design  on  the  basis  of  a 
certain  modification  of  the  O-optimality  criterion.  In  this 
dissertation  we  shall  extend  Khuri's  design  criterion  to  the  cases  where 
we  approximate  the  mean  response  function  given  in  (2.1)  with  spline 
interpolating  polynomials  (with  single  input  variable)  and  Lagrange 

interpolating  polynomial  (with  two  input  variables). 

I* 
Let  R   denote  the  experimental  region  for  the  model 


y  =  f  (x  ;  9)  +  e 


2.32) 


We  suppose  that  R   is  of  the  form 


x,, 


,  xk)  :  a.  <  xi  <  bi  ;  i  =  1  =  2,  ....  k 


By  a  change  of  variables  of  the  form 


2x  .  -a  .  -b  . 

1     1 


'2.33) 
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k  k 

the  region  R   is  changed  to  the  region,  R  =  (t_  :-l  <  t.  <  1,  i  =  1, 

2,  ...,  k}.   If  we  substitute  _t_  for  _x_  in  model  (2.32)  using  the  linear 

transformation  (2.33)  we  obtain  the  model 


y  =  n(t;  8)  +  e. 


2.34) 


Suppose  the  interpolating  polynomial  g(t_;  _9)  is  adequate  for 
approximating  the  response  function  n(t_;  9_)  given  in  (2.34)  over  the 
region  R  and  for  all  Q_  in  some  parameter  space  9;  that  is,  g(t_;  _9_)  is 
approximately  equal  to  n(t_;  _e)  with  some  small  error  of  approximation 
over  the  region  R  and  for  all  _9  in  0.  Thus,  an  approximate 
representation  of  model  (2.34)  is  given  by 


y  =  g(t;  9)  +  e 


12.35) 


A  criterion  for  the  selection  of  a  design  for  estimating  _9  in  model 
(2.35)  which  does  not  depend  on  the  parameter  vector  9  is  given  by  Khuri 
(1982).  This  criterion  is  restricted  to  the  case  of  a  single  input 
variable.  We  shall  extend  Khuri's  criterion  to  the  general  case  where 
we  have  k  input  variables  and  later  in  this  section  we  consider  the 
cases  where  k=2  for  Lagrange  interpolating  polynomial  and  k=l  for  spline 
interpolating   polynomial. 

Let     _tj_,     t_2,      ...,     t^     be     n-points      (n      >     p)     with     coordinates 


,tn,   t12, 


•'     t1k^'     (t?1>     t??>      •••'     t 


21'     L22' 


■Zk),     •••,     (tnl,    tn?,     ...,    tnk 

I, 
over  the  region  R  =  (t_  :  -Kt.<l,  i  =  1,  2,  ...,  k}.   Each  element  ir 
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the  coordinates  corresponds  in  a  one-to-one  manner  through  (2.33)  to 
the  points  of  the  design  D_  =  [x^],  i=l,  2,  ...,  n;  j=i,  ?,  ...,  k  with 
2Lj »  1=1'  2>  •••'  n>  belonging  to  the  experimental  region 

Rx  =  {x:a..<  Xj  <  b . ,  j  =  1,  2,  ...,  k}. 

In  the  design  criterion  by  Khuri  (1982),  some  function  of  the 
determinant  |M_(D_,  ±)  |  is  made  as  large  as  possible  with  respect  to  t,, 
t_2  ....  t^,  where  M_(_D,  _9_)/a2  is  the  information  matrix  of  order  pxp  with 
respect  to  the  model  (2.35)  and  M  has  the  form 


M(D,  9)  =  l     p(t.  ,  9)  r'(t, ,  9) 
1  »1   " n     -  -i  - 


:2.361 


where  the  pxl  matrix   r_(  y ,    Q)    is  the  partial    derivative  of  g(t_i;  _6)   with 
respect    to    9_.      In   matrix    formula,    as   will    be   seen    later,    we   can   write 


Idi  .  D 


U(9)v(^. 


•2.37) 


where  _U(9_)  is  a  pxq  matrix  which  does  not  depend  on  the  design  points 
and  v.(t_.j)  is  a  qxl  matrix  not  involving  _9;  the  value  of  q  will  be 
specified  depending  on  the  form  of  the  interpolating  polynomial.  From 
(2.36)  and  by  using  (2.37),  we  obtain 


M(D,  9)  =  U(9)  I      v(t. )v'(t. )  U'(6) 
1*1    1     1 


'2.33) 


If  we  let 
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A(t_r  t2,  . 


..,  t  ' 

— n  • 


[  v(t.)v'(t. 
i=l   -1  ~  -1 


(2.39] 


then  from  (2.38)  and  (2.39)  we  have 


|M(D,  8)|  =  |U(e)A(tlf  t2,  ....  t, 


)U' 


2.40) 


Let  X  (t.,  t9,  ...,  t  )  and  X^  .  (t,  ,  t0,  ...,  t  )  denote  the  smallest 
(J  — 1  —2      — n        q-1  — 1  —2      — n 

and  largest  eigenvalues  of  A(t_,  t_  ...,  t  ),  respectively.  It  will 
be  shown  later  that  A(t_  x  ...,  t  )  is  nonsingular,  and  hence 
positive  definite,  so  that  the  eigenvalues  are  positive.  From  (2.40)  we 
have 


xP(tr  ...,  tn)|u.(e_)y.'(e_)  |  <  |m(d,  e )  |  <  xp    (t 


q-r-r 


t^)|u(e)u'(e; 


2.41: 


If  the  pxq  matrix  _U(_9_)  is  of  rank  p  then  the  pxp  matrix  U(_9)U '  (_8)  has 
a  nonzero  determinant.  Thus,  we  can  write  (2.41)  in  the  form 


|M(2.  Dl  =  [^(t^ 


-,  t  )  +  (1  -  y)X 


q-1^1' 


tn)J|u(i)U'(_e) 


'2.421 


where  0<y<1.   If  we  assume  that  y  has  a  uniform  distribution  with  the 
densi  ty 


1 


V  Y)  =' 


0<Y<1 

otherwi  se 
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then  from  (2.42),  the  expected  value  of  |M_(D_,  _9)  |  with  respect  to  y   is 
equal  to 


where 


;y[|M(D,  6)|]  =  *(t1,  t2,  ....  ^  )  I  U(_9)U '  (8 )  |       (2.43) 


*(tr  12,    ...,  tj    -   [Xj(tlf  12>    ....  t„)    +  Xj_l(il,  l2,    ....  tfl)]/2 

(2.44) 


The  parameter-free  design  is  the  one  which  maximizes  E  [|M(D,  8)|]  with 
respect  to  t^,  jtg,  ....  t_.  This  can  be  carried  out  by  using  the 
computer  program  SEARCH,  which  was  written  by  Conlon  (1979)  and  is  based 
on  Price's  (1977)  Controlled  Random  Search  procedure,  to  maximize 

*(ii>  i2'  •••'  V  in  (2-44)- 

2.3.1  Designs  Obtained  through  Lagrangian  Interpolation  of  Nonlinear 
Models  in  the  Case  of  Two  Input  Variables 

In  this  section  we  develop  a  method  for  constructing  an  adequate 

approximation  to  the  mean  response  function  given  in  (2.32)  with  a 

Lagrange  interpolating  polynomial.   We  shall  restrict  our  consideration 

of  model  (2.32)  to  the  case  of  two  input  variables. 
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2.3.1.1  Lagrangian  Interpolation  of  the  Mean  Response 

Consider  the  model  (2.32)  for  the  case  of  two  input  variables  (k=2) 

y  =  f(xr  x2;  _9)  +  e   ,  (2.45) 

with  the  experimental  region  R  =  {(x,,  x„)  :  a,<x,<b.,  a?<  x9<b9}.  If 
we  substitute  ti,  to  for  Xi  and  Xo  in  model  (2.45)  using  the  linear 
transformation  given  in  (2.33)  with  k=2  we  obtain  the  model 


y  =  n(t1,  t2;  _8)  +  e   ,  (2.46) 

2 

with  the  region  R  =  {(t,,  t«)  :  -Kt,<l,  -l<t?<l}.  We  assume  that 

(a)   n(t j_ ,  t2",  _9)  has  continuous  partial  derivatives  up  to  order  mi  +  1 

with  respect  to  t^  and  order  m2  +  1  with  respect  to  to  over  the 

2 

region   R       for   all   _9^,   where   m^   and   m2  are   such   that    (m^  +   I)  (mo  + 

1)    is   a  positive  integer  greater  than  or  equal    to  p,   the  number  of 
parameters   in  model    (2.46),   and  are   large  enough  so   that 

(m,+l)  (m2+l) 

max      \r\  (t,,   t.;  _9 )  |  max      |nt  (t,  ,   t9;    8)| 

-Kt^l     ll  i       c  -Kt.<l     Z2  l       l    ~ 

-Kt2<l  -Kt2<l 


m,  m 

(m^l)!      2   l  (m2+l)l      2   c 


(itli+D  (m?+l)  (2'47; 

max      |n.  (t,  ,   t9:    6)n  (t,  ,   t0;   6) | 

-Kt^l       1  l       C    ~    ZZ  l       l    ~ 

-Kt2<l 

+  _ <    5 


(m1+l)!(m2+l)!    2 
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for  all  9  in  some  parameter  space  0,  where 


(m,+l)  m,+l  m.+l 

n*    (t, ,  t,;  e)  =  3   n(t. ,  t,;  e)/3t, l 


't,    v"r  -2' 


1'  '2'  ^"""l 


(m«+l)  m?+l  m„+l 

nt    (tj,  t2;  _e)  =  3    n(tls  t?;  _9)/3t2 


and  6  is   a  given  positive  constant. 
(b)        For    any    set    of    (irij+1)  (mg+l)    distinct    points    (yQ,    Zg),    (yg,    Zi), 


•••,     (yn.    0>     (yi>    zn),     ....     (y_   .    z_   ),    such    that 


■i<y0< 


and 


•l<zn<z1 <. . .<z     <1     where   mi    and   m9   are   the 
u     l  nip  J-' 

integers   defined   in  Assumption    (a),   the  px   [(m^+1) (mo+1)]  matrix, 


y.<...<y  <1 

y  1    'in. 


U(£)  =  [U(y0,  V  9),  U(y0,  zr  8),  ...,  U(ym  ,  z^  ,  9)]   ,  (2.48) 


has  rank  p,  where  U_(y.  ,  z.,  _9_)  is  the  partial  derivative  of 

My...  z.;  8_)    with    respect  to    _9  (i  =0,  1 ,  . . . ,  m  •  j  =0, 

1,  ...,  m2). 

We    consider    the   points    with  coordinates    (yg,    Zg),    (y«,    z-,), 

...,    (y     ,    zm   )    described    in    Assumption    (b)    with   y-    and    z,   being 
12  J 

of  the  Chebyshev  type  given  in  (2.14)  with  a^  =  -1,  bi  =  1,  a?  = 
-1  and  b2  =  1.  Substituting  t^,  t2  for  x^  and  X2  in  (2.10)  using 
the  linear  transformation  given  in  (2.33)  with  k=2,  the  Lagrange 
interpolation   polynomial    becomes 


where 


40 


g(tr  t2; 


I        n(y:-,  z.:e_)}..(tr  t?\ 
0<i<m, 


0<j<m 


2.49) 


VV  V 


=    n   — — -   n   -=-— - 

k=0  yi~yk  u=0  zj"zu 


k*i 


u*j 


i=0,   1. 
j=0,   1, 


'2.50) 


The  error  resulting  from  the  use  of  Lagrangian  interpolation  is,  by 
definition,  the  difference  between  n(t.,  t?;  9)  and  g(t.  ,  t  •  a),  the 
upper  bound  of  this  error  can  be  obtained   from   (2.15)   and  is   given  by 


max      | e ( t . ,  t „ ) |   < 
■Kt^l  c 

•Kt2<l 


(m.+l) 

max      \r\  (t.  ,   t    : 

•Kt,<l       1  l       d 

•l<t2<l 


^+1)12 


(m2+l) 
max        |    r\  (t:  ,   t0; 

■l<t.<l         z2 

■Kt2<l 


I-   <-2'  -' 


(m2+l)!      2 
(m1+l) 


(m2+D 


™x      |n  (t,,   t?;    8)tv      "  (t    ,   t9. 


■Kt?<l 


(m1+l)!(m2+l)!2 


m.  -Hrip 


'2.5i: 
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2.3.1.2  Designs  Based  on  Lagrangian  Interpolation 

Consider  the  Lagrange  interpolating  polynomial  g(t^,  to;  _8)  given 
in  (2.49)  with  interpolation  points  being  of  the  Chebyshev  type  given  in 
(2.14)  with  a^  =  -1,  b^  =  1,  ao  =  -1  and  b2  =  1-  The  positive  integers 
mi  and  m2  are  determined  such  that  g(ti,  to;  _9_)  is  approximately  equal 

to  n(t.,  t?;  _8_)   with  an  error  of  approximation  not  exceeding  5  over 

2 

the  region   R   and  for  all  _9_  in  0.   An  approximate  representation  of 


model  (2.46)  is  given  by 


y  =  g(tlf  t2;  _e)  +  e 


'2.52) 


Let  (tji,  t^),  (t21'  t22^'  "•'  ^nl'  tn2^  be  the  P01'nts  in  tne 
p 
region  R.  =  {(t,,  to) : -l<t,<l,  -l<tp<l}  which  have  a  one-to-one  corre- 
spondence through  (2.33)  with  the  points  D_  =  [(xii,  x12) »  (*21'  x22^' 

p 
...,  (xnl,  xn2)]'  in  the  region  R~  =  {(x^,  x2)  ta^x^b,,  a2<Xp<b2}.   We 

consider  Kt^  ,  e_)  in  (2.36)  since  it  is  the  partial  derivative  of  g(til, 
ti2;  D  with  resPect  to  A'  therefore,  from  (2.49)  we  can  write  £_( t_j ,  _9) 
as 


3ri(y  ,  z.  ; 
r(t,,  9)  -I  U»      J 


0<u<m, 
0<j<m, 


39 


luj(tn,ti2)  ,  1-1.  2. 


'2.53) 


From  (2.37)  and  (2.53),  the  px[(m1  +  l)  (m2  +  l)]  matrix  U_(_9_)  can  be  written 
as 
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U(9)    = 


3n(y0,   2Q;   9)       3n(yQ,   z^   9) 


3ri(yn 


z      ;    9N 
m„     — ' 


36 


39 


2.54) 


and    is    of    rank    p    by    Assumption    (b)    and  _v(_t_n- )    is    an    [(fth  +  I)  (riip+l)]   x    1 
matrix  of  the   form 


v(t-)   =   [lnn(t.,,   t.0),   lm(t.1s   t.0),    ...,   1  ft.,,   t.    ' 

— v— i  ;        L   OCT   il'     i2y       01v   il'     i2  m,m?v   il'     i2 


2.55) 


If  n  >  (m,+l)(m„+l)  and  at  least  (m^+1)  of  t^,  ...,  t^.are  distinct, 

^2,   •••.  t_2   are  distinct,   then  the 


and  at  least  (m0+l)  of  t 


t  )  of  the  form 
— n 


[(m,+l)(m2+l)]  x  [(m,+l)(m2+l)  ]  matrix  A(t,,  t_?, 
(2.39),  where  _v ( t- )  is  given  in  (2.55),  should  be  of  full  rank  and  thus 
nonsingular.  To  show  this,  suppose  that  A(t_,  ,  t~,  ...,  t  )  is  of  rank 
less  than  (m^+1) (n^+l)  tnen  tne  (m^+lMiT^+l)  x  n  matrix 
A  =  [^.(IJ'  li^jp) '  •••>  lit  )]  is  of  the  same  rank  as  A(t_,,  _t?,  ..., 
t^].  In  this  case  there  must  exist  a  set  of  a  •,  u=0,  1,  ...,  m,; 
j  =0,  1,  ...,  mo,  not  all  equal  to  zero  such  that 


I         I     a    .1    .(t.,  ,  t.0)   =  0       ,   i=l,   2, 
L         L       u/j   ui v   l 1'     i2'  '  '      ' 


,   n 


u=0     j=0 


(2.56) 


on   R' 


(t,,   t2):    -Kt^l,   -l<t2<l)}.     Since  we  can  write   l^-U^,    ti2' 

m.       m„ 

lu(tn)lj(ti2)     so     (2.56)     becomes       [         [  a.l(t.  .)1  .  (t.  -)   «  0, 

u=0     j=0  J                 J 
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1=1,  2, 


Let   Cj(tn: 


"1 
I  a 
u=0 


.1  t., 

UJ  U v  1  1 


It  follows  that 


But 


I     C(t.,)1  .(t.0)  =  0  for  each  t, ,  in  [-1,1],  1-1,  ?,  ,.,,  n. 

the  Ij's  are  linearly  independent  on  the  interval  [-1,1]  (see  'Khun' 

(1932)).  Thus,  Cj(ti  j_)  =  0  for  each  j=0,  1,  ....  m2,  and  each  tn,  1=1, 

2,  ...,  n,  in  [-1,1].   But  then  for  each  j,  j=0,  1,  2,  ...,  m?, 

ml 

y  a  .1  (t.,)  =  0,  1=1,  2.  ...,  n,  and  by  the  linear  independence  of 

the  lu's,  au  •  =  0  for  all  u=0,  1,  ...,  m^  and  j=0,  1,  ...,  m?,  which  is 
impossible  since  we  assume  that  aui-'s  (u=0,  1,  ...,  m,;  j=0,  1,  ...,  m2) 
are  not  all  equal  to  zero. 

We  now  have  U_(_9_)  a  matrix  of  rank  p,  and  thus  U_( 6_)U_'  ( 9_)  has  a  non- 
zero determinant.  We  have  also  shown  that  the  matrix  Aft,,  t0,  ...,  t  ) 

— 1  —  l  — n 

is  of  rank  (mi+1) (trio+l) ,  3nd  thus  nonsingular.  From  Section  2.3,  the 
parameter-free  design  maximizes  ip(t_,  ,  t_?,  ...,  t  )  with  respect  to  ti, 


t^>,  . ..,  t^  where  ty(t±,   j^?, 


is  of  the  form 


*(t,,  ...,  t  ' 
— 1 '     -n 


■  Wh* 


t 

-n 


+  \> 


(m1+l)(m2+l)-l 


-1* 


t 


J/2 

'2.57: 


where   xn(li  >  •  ■ 


t 

-HI 


and   X  (t_  , 

(m1+l)(m2+l)-l 


t  )  denote  the 
-n ' 


t 

^n 


smallest  and  largest  eigenvalues  of  A(t_,  ,  t_?, 

We  have  provided  the  methodology  for  obtaining  the  optimal  design 
based  on  Lagrange  interpolating  polynomial  for  the  case  of  two  input 
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variables.  It  was  shown  that  the  number  of  design  points  is  at  least 
equal  to  (m^+lH^+l)  •  This  number  tends  to  grow  rapidly  as  mi  or  m? 
become  large.  In  this  case,  it  becomes  vary  difficult,  and  practically 
impossible,  to  determine  optimal  design  settings  using  Price's  computer 
algorithm,  or  any  other  similar  algorithm,  due  to  the  high 
dimensionality  of  the  optimization  problem. 

2.3.2  Designs  Obtained  through  Spline  Interpolation  of  Nonlinear  Models 
In  this  section  we  develop  methods  for  constructing  designs 
assuming  the  true  mean  response  given  in  (2.32)  can  be  adequately 
approximated  with  a  linear,  a  quadratic  and  a  cubic  spline  in  the  case 
of  a  single  input  variable.  Each  approximation  will  be  used  to  obtain  a 
design  for  parameter  estimation  which  does  not  depend  on  the  parameter 
vector  9_. 

2.3.2.1  Spline  Interpolation  of  the  Mean  Response 

Consider  the  model  (2.32)  for  the  case  of  a  single  input  variable 
(k-1) 

y  -  f(x;  _6)  +  e  (2.58) 

with  the  experimental  region  R  =  {x:  a<x<b}.  If  we  substitute  t  for  x 

in  model  (2.58)  using  the  linear  transformation  given  in  (2.33)  with  k=l 
we  obtain  the  model 
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y  =  n(t ;  9)  +  e 


2.59] 


the  region  R   is  reduced  to  the  interval   -l<t<l, 
3     x 


Consider  the  knots  tq,  t,,  ...,  t^ 


where  -1  =  xn  <  x      <  . . .  < 


1 


Th  <  Th+1  a  !•  A  spline  of  degree  m  (with  m-1  continuous  derivatives  at 
the  knots  Tj,  i=l,  2,  ...,  h)  which  approximates  n(t;  9_)  can  be 
obtained  by  substituting  t  for  x  in  (2.17),  we  obtain 


9(t;  9)  =  I     a.t1  +  I    b  (t  -  t.j;  , 
1*0  n     i=l 


2.60) 


where 


(t-T, 


,   t>T. 


It  -r.K 


otherwi  se 


,  j-1,  2,  ...  (2.61! 


The  coefficients  ag,  a^,  ...,  a  ,  bi,  b?,  ...,  b^  in  model  (2.60)  can  be 
obtained  by  solving  the  interpolation  problems  in  (2.21),  (2.25)  and 
(2.29)  for  linear,  quadratic  and  cubic  splines,  respectively.  We  note 
that  the  coefficients  a^  ,  b-,  i=0,  1,  2,  ...,  m;  j=l,  2,  ...,  h  are 
functions  of  _9_. 

Consider  the  model  in  (2.59).  We  assume  that 

(a)   n(t;  9_)  has  continuous  derivatives  up  to  order  m+1  with  respect  to 

t  over  the  interval   -l<t<l   for  all  _9_,  where  m  is  the  degree  of 

the  spline  function  and  is  such  that  m+h  +  1  is  a  positive  integer 

greater  than  or  equal  to  p,  the  number  of  parameters  in  model 
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(2.58)  and  h  is  the  number  of  interior  knots  and  is  large  enough 
so  that 

2  (2) 

(At)       max    |nv    '  (t ;  _9 )  |  /   8  <   5,   for   linear  spline,    i.e.,   m=l  , 

-Kt<l 

(2.62) 

or 

3  (3) 

(At)       max    |nv    ' (t  ;   6)|/8   <   6,    for  quadratic   spline,    i.e.,   m=2, 

-Kt<l 

(2.63) 

or 

5(At)       max    |n^(t:  _9)|/384  <   <5,    for  cubic   spline,    i.e.,  m=3, 
-Kt<l 

(2.64) 

for  all  _9_  in  some  parameter  space  9,  where  At  =  t.  -  t.  , 
i  =  0,  1,  — ,  h,  is  the  common  distance  between  the  knots, 
n(j)(t;  9)  =  3jn(t;  _9)/3tJ',  j  =  2,  3,  4,  and  5  is  a  given 
positive  constant.  We  need  the  quantity  m+h+1  to  be  a  positive 
integer  greater  than  or  equal  to  p  because  of  the  parameter  esti- 
mation purpose;  since  m+h+1  is  the  number  of  a's  and  b's  in  (2.60) 

and  these  a's  and  b's  are  functions  of  9  =  (9,,  ...,  9  )'.   if 

1       p' 

m+h+1  <  p  we  cannot  uniquely  estimate  9   9  .  The  assumption 

that   n(t:  8_)   has  continuous  derivatives  up  to  order  m+1  with 
respect  to  t  is  needed  for  the  purpose  of  obtaining  the  upper 
bound  of  the  error  of  interpolation. 
The  px  (m+h+1)  matrix 

U(6)  =  [1^(8),  Ujte),  ...,  J4„+h(8)]  (2.65) 


47 


is   of   rank   p,   where   U^  (e_)   =  3a./36_,  i  =  0,  1,  ...,  m, 

U  .(6)  =  9b. /96   i  =  1,  2,  ...,  h  and  a,,  b.-,  i=0,  1,  ...,  m, 
— m+j  —     J  -.  1   J 

j  =  1,  2,  ...,  h,  are  the  coefficients  in  model  (2.60). 

The  error  resulting  from  the  use  of  spline  interpolation 
polynomial  is  the  difference  between  n(t;  9_)  and  g(t;  9_) .  The 
upper  bound  of  the  errors  can  be  obtained  by  substituting  t  for  x 
in  (2.23),  (2.27)  and  (2.31)  using  the  linear  transformation  given 
in  (2.33)  with  k*l  resulting  in 


max  |g(t; 
•l<t<l 


n(t;  e_)|  <  (Ax)2  max  |n^(t;  9_)|/3, 
-Kt<l 


for  linear  spline, 


(2.66) 


max  |g(t ;  6_)  -  n(t ; 
•Kt<l 


<  (At)3  max  h(3)(t 
-Kt<l 


/8, 


'2.67) 


max  |g(t ;  £] 
■Kt<l 


for  quadratic  spline, 


n(t;  9_)  |  <  5(Ax)4  max  |n^(t 
-Kt<l 


9)1/384, 


'2.63) 


for  cubic  spline 


where  Ax  and  n'J^(t;  _8 ) ,  j  =  2,  3,  4,  are  defined  in  Assumption 
(a). 


2.3.2.2  Designs  Based  on  Spline  Interpolation 

Consider  the  spline  interpolation  polynomial  g(t;  _9_)  given  in 
(2.60)  with  h  interior  knots.  The  positive  integer  h  is  determined  by 
requiring  the  upper  bounds  in  (2.66)  or  (2.67)  or  (2.63),  for  linear, 
quadratic  and  cubic  spline,  to  be  less  than  some  specified  small  value 
5  >  0  for  all  •  9_  in  some  parameter  space  0.  Thus,  an  approximate 
representation  of  model  (2.59)  is  given  by 


y  =  g(t;  9)  +  e 


'2.69) 


Let  tj_,  tp,  ...,  tn  be  the  points  on  the  interval  -l<t<l  which 
correspond  in  a  one-to-one  manner  through  (2.33)  to  the  points  of  the 
design  D_  =  [xp  X£,  ...,  xp]'  with  xi  (i=l,  2,  ...,  n)  belonging  to  the 
experimental  region  R  =  {x:  a<x<b}.  From  (2.60)  we  can  write  _r(t.,  j3) , 
where  r_(t^ ,  9_)  is  the  partial  derivative  of  g(t;  _e_)  with  respect  to  _9_, 
as 


^i'i)  -  j0  Wu^l  +  u_-l  W-i 


uJ  +  ' 


i  =  1,  2, 


From  (2.37)  and  (2.70)  we  obtain  U_(_9)  and  v_(t.-)  of  the  form 


,  n 


'2.70) 


U(9) 


3     3 


3_     3L  3. 

'  Ta"l. 


■39^0'  39  1'  ••''  39m'  39  T  '**'  IF  h 


:2.7i: 
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v(t.)  =  [1,  t.,  t^  . 


^  -  v+> 


2.1V, 


the  matrix  U(8)  is  of  order  px  (m+h+1)  and  of  rank  p  by  Assumption  (b' 


and  v_(tj)  is  an  (m+h+1  )xl  matrix. 


Now,  if  n>m+h+l  and  at  least  m+h+1  of  the  design  points  are 
distinct  (this  is  needed  so  that  the  a's  and  b's  in  (2.60)  are  uniquely 
estimated),  then  the  (m+h+1)  x  (m+h+1)  matrix 


A(tr  t2,  ...,  tn)  =  I     v(t.)v'(t.)    , 


'2.731 


should  be  of  rank  m+h+1  and  thus  nonsingular.   To  show  this,  suppose 
.,  t  )  is  of  rank  less  than  m+h+1,  then  the  (m+h+l)xn 


that  A(tp  t2, 
matrix  G_  =  [v_(ti), 


v_(tp)]  is  of  the  same  rank  as   A(t.,  t„, 


Therefore,  there  must  exist  a  set  of  a.-,  i  -   0,  1,  ...,  m+h , 


not  all  equal  to  zero  such  that 


Tl 

m+h 

I 

a 

t! 

+     )' 

i«0 

j 

i  =m+l 

N(tj 


1  -m'  + 


,3=1,  Z.. 


2.74) 


m     .    m+h 
From  (2.74)  we  conclude  that   I     a.t1  +   £   a.(t  -  t.   )m,   which  is 

i  =0       i  =m+l 
a  polynomial  of  degree  m,  has  n  roots,  namely,  ti,  to,  ...,  t_,  which  is 
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impossible  since  n>m+h  +  l>m  and  at  least  m+h+1  of  ti,  to,  ..., 

tn  are  distinct  (a  polynomial  of  degree  m  has  at  most  m   distinct  roots). 

We   have   thus   shown   that   the   (m+h+1)  x  (m+h+1)   matrix 

A(t.,  t   ...,  t  )   is  of  rank  m+h+1.   Therefore,  it  is  positive 
1   2      n  . 

definite  and  by  Assumption  (b),  JJ_(_9)  is  a  matrix  of  rank  p,  hence  the 
pxp  matrix  U_(9_)U_' (9_)  has  a  nonzero  determinant.  Consequently,  from 
Section  2.3,  the  parameter-free  design  is  the  one  which  maximizes 
<P(t,,  t  ...,  t  )  with  respect  to  t^,  t?J  ...,  tn  and 
i>[t,,   t?,  ...,  t  )  is  of  the  form 

♦  (tlf  t2,  ....  tn)  -  [xP(tl.  t2,  ...,  tn)  +  XP+h(tr  t2.  ...,  tn)]/2 

(2.75) 
where  ^(t^  t^,  ...,  t  )  and  *  ^(t,,  U,  "•'  tn^  denote  the  smallest 
and  largest  eigenvalues  of  A(t, ,  t„,  ...,  t  ). 

We  note  that  in  order  to  be  able  to  generate  a  satisfactory  distri- 
bution of  information  throughout  the  interval  [-1,1],  the  optimal  design 
obtained  on  the  basis  of  spline  interpolation  is  subject  to  the  con- 
straint that  at  least  one  design  point  is  allocated  in  each  subinterval. 


2.4  Exampl es 
We  now  present  two  examples  for  obtaining  parameter-free  designs 
based  on  Lagrangian  interpolation  (for  the  case  of  two  input  variables) 
and  spline  interpolation  polynomials  of  degrees  1,  2  and  3  (for  the  case 
of  a  single  input  variable). 
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Example  2.1  Consider  the  model  used  by  Box  and  Hunter  (1965)  for  a 
certain  chemical  reaction  rate  which  can  be  described  by  the  nonlinear 
model 

f(x1P  x2;  8)  =  9^^/(1  +  8^  +  92x2)   ,        (2.76) 

where  f(x1?  x2;  _9)  represents  the  reaction  rate,  Xi  and  x2  are  partial 
pressures  of  reactant  and  product  respectively,  6^  and  92  are  adsorption 
equilibrium  constants  for  reactant  and  product  respectively,  and  83  is 
the  effective  reaction  rate  constant.  Suppose  that  a  prior  knowledge  of 
the  parameters  9^,  92  and  83  is  available  in  the  form  of  known  lower  and 
upper  bounds.  Thus,  Vj  <91<w1>  v2<92<w2  and  v3<93<w.  for  some  known 
values  vi  and  w^i  =  1,  2,  3)  and  suppose  that  the  experimental  region 
is  Rx~  =  {(Xp  x2):  0<x1<bp  n<x2<b2}  where  b1  and  b2  are  some 
positive  values  specified  by  the  experimenter  (in  this  example  we  choose 
bl  =  3  and  b2=2).  Using  the  linear  transformation  given  in  (2.33)  with 
k=2,  a^O  and  a2=0  model  (2.76)  becomes 

n(t1,  t2;  9)  =  e1e3b1(t1+i)/[2  +  e1b1(t1+l)  +  92b2(t2+i)]   , 

(2.77) 

2  ? 

the  region  Rx  is  transformed  to  R^  =  {(t^  tg):  -Kt,<l,  -Kt2<l}.  The 

Lagrangian  interpolation  accuracy  is  determined  by  using  the  equation 

given  in  (2.51).  It  can  be  shown  (see  Appendix  1)  that  for  all  8  in  0, 

the  upper  bound  for  the  error  is  given  by 
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m?+l     nu+1 
,  w.w,,b1(t1+l)w?         by 

max     |e(t,  ,t0)  |   <  -1-      max     I — =— i-= — = =•! 

Mt   <1  i      '  [TIo      ,    .       -■  m,,  +  <i 

:{<?i<;  2   :l<ti<i  [2+v1b1(t1+D+v2b2(t2+i)] 2 


m,+l       m.+l 
x                        w1         w3b1       [2+w2b2(t2+l)] 
max     I -s-l 

2       -i^l    [2+v1b1(t1  +  l)+v2b2(t2+l)]    l 


— rr-       max 
2ml4fll2  -Kt^l 

-Kt2<l 


m,+2  2  m, +2  m?+l     m?+l 

wl      w3bl       (ti+1)w2        b2       [2+w2b2(t2+l)] 


;2+v1b1(t1+l)+v2b2(t2+l) 


m,+nip+4 


2.78) 


Suppose  now  we  require  that  the  upper  bound  of  the  error  in  (2.73) 
(with  bj=3  and  b2=2)  be  less  than  or  equal  to  5  =  .05.  We  note  here 
that  this  upper  bound  involves  the  upper  and  lower  limits  of  8-,,  82  and 
83.  Suppose  that  the  lower  and  upper  limits  of  9,,  8?  and  93  are  such 
that  .5  <  61  <  1,  .7  <  92  <  2  and  .4  <  93  <  .8.  These  values  are 
obtained  from  the  approximate  95%  confidence  region  of  _9  which  was 
constructed  using  ten  response  values  that  were  generated  according 
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to  the  model  in  (2.77)  with  9^  82,  63  and  the  error  variance,  a2,  taken 
as  .6,  1.2,  .5  and  .0006,  respectively.  The  ten  design  points  were 
chosen  from  the  region 


R^  =  {(tl,   t2):  -1  <  tj_  <  1  ,  -1  <  t2  <  I] 


To  compute  the  upper  bound  of  the  error  we  fix  m^  and  m2  and  carry  out 
the  maximization  as  in  (2.78)  over  the  region  R^  =  {(t,,  t?): 
-1  <  t.  <  1,  -1  <  t-  <  l}.  This  maximization  is  accomplished  by  the 
use  of  the  computer  program  SEARCH,  which  was  written  by  Con!  on 
(1979).  By  computing  this  upper  bound  for  several  assignments  of  mi  and 
2  we  obtain  the  values  of  m^  and  m2  to  satisfy  the  <5  =  .05 
requirement.  These  values  are  mj_  =  12  and  m2=4  and  the  value  of  the  upper 
bound  is  .04816541.  The  Chebyshev  points  corresponding  to  the  values  of 
m1  and  m2  are  obtained  from  (2.14)  with  a-p-1,  a2=-l,  b,=l  and  b2=l. 
These  points  are  given  in  Table  1.  We  use  these  Chebyshev  points  to 
form  the  coordinates  of  points  in  two  dimensions,  we  obtain  (mi+l)(m2+l) 
=  55  points.  These  coordinates  are  given  in  Table  2.  We  choose  n,  the 
number  of  design  points,  to  be  equal  to  65  (here  all  sixty-five  points 
are  distinct).  The  vector  _u ( t^ )  as  well  as  the  matrix  A(tt,  to,  ..., 
t_£5)  can  be  obtained  from  formulas  (2.55)  and  (2.39),  respectively, 
which  leads  us  to  the  maximization  of  <|>(ti,  t2,  •••,  tec)  in  (2.57)  with 
p=3.  Since  the  number  of  design  points  required  for  this  example  is 
computationally  intractable,  the  actual  determination  of  the  65-pt. 
design  was  not  attempted. 


m 
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Table  1 
Chebyshev  Points  for  Example  2.1 


Based  on  m^  =  12  Based  on  m?  =  4 


.99271 

.95106 

.93502 

.58779 

.82298 

0 

.66312 

-.58779 

.46472 

-.95106 

.23932 

0 

-.23932 

-.46472 

-.66312 

-.82298 

-.93502 

-.99271 

55 


Table  2 


Coordinates  of  Points   in   Two-Dimension 
Obtained   from  Chebyshev  Points   for  Example  2.1 


(.99271, 

.95106) 

(    .46472, 

0      ) 

(-.46472, 

-.95106) 

(.99271, 

.58779) 

(    .46472, 

-.58779) 

(-.66312, 

.95106) 

(.99271, 

0     ) 

(    .46472, 

-.95106) 

(-.66312, 

.58779) 

(.99271, 

-.58779) 

(    .23932, 

.95106) 

'-.66312, 

0     ) 

(.99271, 

-.95106) 

(    .23932, 

.58779) 

(-.66312, 

-.58779) 

(.93502, 

.95106) 

'    .23932, 

0     ) 

-.66312, 

-.95106) 

(.93502, 

.58779) 

.23932, 

-.58779) 

-.82298, 

.95106) 

(.93502, 

0     ) 

.23932, 

-.95106) 

-.82298, 

.58779) 

(.93502, 

-.58779) 

0     , 

.95106) 

-.82298, 

0     ) 

(.93502, 

-.95106) 

0     , 

.58779) 

-.82298, 

-.58779) 

(.82298, 

.95106) 

0     , 

0     ) 

-.82298, 

-.95106) 

(.82298, 

.58779) 

0     , 

-.58779)            ( 

-.93502, 

.95106) 

(.82298, 

0     )           ( 

0      , 

-.95106)           ( 

-.93502, 

.58779) 

(.82298, 

-.58779)            ( 

-.23932, 

.95106)           ( 

-.93502, 

0     ) 

(.82298, 

-.95106)            ( 

-.23932, 

.58779)           ( 

-.93502, 

-.58779) 

(.66312, 

.95106)            ( 

-.23932, 

0     )           ( 

-.93502, 

-.95106) 

(.66312, 

.58779)            ( 

-.23932, 

-.58779)            ( 

-.99271, 

.95106) 

(.66312, 

0     )            ( 

-.23932, 

-.95106)            ( 

-.99271, 

.58779) 

(.66312, 

-.58779)            ( 

-.46472, 

.95106)            ( 

-.99271, 

0     ) 

(.66312, 

-.95106)            ( 

-.46472, 

.58779)            ( 

-.99271, 

-.58779) 

(.46472, 

.95106)            ( 

-.46472, 

0      )            ( 

-.99271, 

-.95106) 

(.46472, 

.58779)            ( 

-.46472, 

-.58779) 
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Example  2.2  In  this  second  example  the  parameter-free  designs  are 
obtained  through  linear,  quadratic  and  cubic  spline  interpolating 
polynomials  used  for  approximating  the  nonlinear  model.  We  consider  the 
model  used  by  Draper  and  Smith  (1981).  The  model  is  taken  from  an 
investigation  performed  at  Proctor  and  Gamble.  The  investigation 
involved  the  relationship  between  Available  Chlorine,  y,  of  product  A 
and  the  length  of  time  of  manufacture,  x.  It  is  required  that  product  A 
must  have  a  fraction  0.50  of  Available  Chlorine  at  the  time  of 
manufacture  and  the  fraction  of  Available  Chlorine  in  the  product 
decreases  with  time.  The  form  of  the  model  is 

y  =  a  +  (.49  -  a)exp(-B(x  -  8) )  +  e   .         (2.79) 

We  suppose  that  the  experimental  region  R  =  {x:  a<x<b}  where  a  and  b 
are  some  positive  values  specified  by  the  experimenter  (in  this  example 
a  =  10  and  b  =  40),  and  let  the  parameter  space  0,  be  such  that 
0.36<a<0.41  and  0.06<3<0.16  (these  values  are  obtained  from  the  95% 
confidence  region  of  a  and  6  in  the  example  by  Draper  and  Smith  (1981, 
p.  486)).  Using  the  linear  transformation  given  in  (2.33)  with  k=l,  a=10 
and  b=40,  model  (2.79)  becomes 

y  =  a  +  (.49  -  a)exp(-8(15t  +  17))  +  e   ,         (2.80) 

where  -l<t<l.  It  can  be  shown  (see  Appendix  2)  that  for  all  8  in  9, 


max  |g(t,i 
■Kt<l 


n(t,i 


<  < 
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^At)~( .0679675)     for  linear  spline 


>t)j(. 1631221197)  for  quadratic  spline 


;at)h(. 0407805)     for  cubic  spline 


(2.81) 
Suppose  now  we  require  that  the  upper  bounds  in  (2.81)  be  less  than  or 
equal  to  <5  -  .05,  by  taking  into  consideration  that  we  have  equidistant 
knots.  Then  the  values  of  At  =  t.-+1  -  ti ,  i=0,  1,  ...,  h,  that  satisfy 
the  <5  =  .05  requirement  are 


'.4  with  maximum  error  =  .0108748  for  linear  spline 
At  =  -^  .5  with  maximum  error  =  .02039   for  quadratic  spline 
^1  with  maximum  error  =  .0407805  for  cubic  spline. 


2.82) 


From  (2.82),  the  number  of  interior  knots  on  the  interval  -Kt<l  for 
linear,  quadratic  and  cubic  splines  are  h=4,  3  and  1,  respectively.  The 
knot  positions  for  linear  spline  are  -1,  -.6,  -.2,  .2,  .6,  1,  for 
quadratic  spline  are  -1,  -.5,  0,  .5,  1,  and  for  cubic  spline  are  -1, 

0,  1.  By  choosing  n,  the  number  of  design  points  to  be  equal  to  m  +  h  + 

1,  where  m  is  the  degree  of  the  spline  functions  (in  our  case  m=l,  2, 
and  3)  and  h  is  the  number  of  interior  knots,  we  obtain  n=6,  6  and  5  for 
linear,  quadratic  and  cubic  splines,  respectively.  The  vector  v(t)  as 
well  as  the  matrix  A(t^,  t?>  ...,  tn)  can  be  obtained  from  formulas 
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(2.72)  and  (2.73),  respectively.  The  maximum  values  of  \|»(ti,  t?,  ..., 
tn)  in  (2.75)  with  p=2  for  linear,  quadratic  and  cubic  splines  are 
178.605,  219.227  and  114.813  and  the  corresponding  optimal  values  of 
t -j  ( i  =  1 ,  2,  ...,  n)  along  with  the  join  points  (interior  knots)  t-  ,  i=l, 
2,  ...,  h,  are  given  in  Table  3. 

For  comparison  purposes  and  to  check  the  efficiency  of  our  optimal 
design,  by  using  _9  =  (a,  0)'  =  (.39,  .10)',  values  of  the  determinant 
of  the  matrix  M_(D_,  9_)  given  in  (2.38)  were  computed  at  both  optimal 
design  and  the  Box-Lucas  design  (see  Appendix  3)  and  are  given  in 
Table  4.  The  efficiency  is  defined  by:   (Lucas  (1976)) 

WMWimal  0esiqnl/nl  ^ 
Efficiency  =  [ 


|M(D,9)D       r/nT 
'—v-- Box-Lucas  '  2 


where  n^,  n2>  (in  this  example  rii  =  6,  6,  5  for  linear,  quadratic  and 
cubic  splines  and  n^  =  2),  are  the  numbers  of  design  points  for  the 
optimal  and  Box-Lucas  designs,  respectively,  and  p  is  the  number  of 
parameters  in  the  model  (in  our  case  p=2).  The  efficiencies  are  also 
given  in  Table  4. 

From  Table  4  we  can  conclude  that  with  respect  to  model  (2.79),  our 
optimal  designs  based  on  the  linear  and  the  quadratic  splines  are  1.357 
and  1.042  times  more  efficient  than  the  Box-Lucas  design.  But  the 
optimal  design  based  on  cubic  spline  has  a  lower  efficiency  than  the 
Box-Lucas  design. 

For  the  purpose  of  comparison,  the  plots  of  the  linear,  quadratic 
and  cubic  splines  are   produced  in  Figures  2,  3  and  4.  For  the  true  mean 
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response  curve  in  Figure  1,  we  set  a  =  .39  and  8  =  .10  (tnese  two  values 
are  the  least  squares  estimates  of  a  and  B  obtained  by  Draper  and  Smith 
(1981)).  The  equations  for  the  linear,  quadratic  and  cubic  splines  that 
provide  the  approximations  to  the  model  in  (2.80)  with  the  maximum 
errors   given   in   (2.82)   are  of  the   forms 

linear   spline: 

g(t;6_)    =   .379523-. 0923504t+.0416674(t+.6)++.0223676(t+.2)+ 
+.01255(t-.2)++.00688758(t-.6)+   , 

quadratic  spline: 

g(t;j9)    =   .416754+. 00793627t+.0630552t2-.0340915(t+.5)| 
-.0151525(t-0)f-. 00724385ft-. 5) J, 

cubic   spline: 

g(t;8_)    =   .408268+. 0261167t+.0157712t2 
-.0217168t3+.0178701(t-0)3    . 

(2.83) 

These  splines  are  constructed  under  the  assumption  that  each  of  them  has 
m-1  continuous  derivatives  at  the  knots,  where  m=l,  2,  and  3, 
respectively.  The  mean  response  values  obtained  from  the  model  in 
(2.80)  (with  a  =  .39  and  0  =  .10),  and  values  obtained  from  linear, 
quadratic  and  cubic  splines  in  (2.83)  are  given  in  Table  5  for  certain 
values  of  t  in  the  interval  [-1,1],  the  plots  of  the  response  values 
(for  each  model)  against  the  t  values  are  given  in  Figures  1,  2,  3,  and 
4,  respectively. 
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Table  3 
Optimal  Designs  and  Interior  Knots.  Example  2.2 


Linear  Spl ine  with 
Interior  Knots  at 
-.6,  -.2,  .2,  .6 


Optimal  Designs 

Quadratic  Spline  with 

Interior  Knots  at 

-.5,  0,  .5 


Cubic  Spline  with 

Interior  Knots  at 

0 


-.6960602 
-.200152 
.1999294 
.5999959 
.9997785 
.9999233 


-.9823732 

-.9458824 

-.0002855 

.4995994 

.9999732 

.9999965 


.  -1 

-.0007636 

.9975316 

.9990895 

1 


Table  4 

Values  of  |_M(£,_9 )  |  from  Optimal  and  Box-Lucas  Designs 
for  Q_  =   (a,    &)'    =    (.39,  .10)' 
and  the  Efficiencies  of  the  Optimal  Design.  Example  2.2, 


|M(D,0)   «..,.. 

'—  —  —  optimal  design 


l^'l)  Box-Lucas'   Efficiency  of 
Optimal  Design 


Linear 
Spl ine 

Quadratic 
Spl ine 

Cubic 
Spl  i  ne 


.387684347 
,247975528 
,169248231 


.069146279 


.076197813 


.085737989 


1.367 


1.042 


.889 
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Table  5 

Mean 

Response  Values 

with  Corresponding 
Mean  Response  Val 

Values  of  t. 
ues 

Example  2.2 

t 

True  Model 

Linear  Spline   Quadratic  Spline 

Cubic  Spline 

-iLQC 

.471873 

.471873 

.471873 

.471873 

-.9 

.460469 

.462638 

.460686 

.460379 

-.8 

.450653 

.453403 

.450760 

.450374 

-.7 

.442205 

.444168 

.442096 

.441726 

-.5L 

.434933 

.434933 

.434692 

.434306 

-.51 

.428674 

.429865 

.428674 

.427984 

-.4 

.423287 

.424797 

.423327 

.422628 

-.3 

.418650 

.419728 

.418684 

.418109 

-.2L 

.414660 

.414660 

.414621 

.414296 

-.1 

.411255 

.411878 

.411136 

.411059 

0QC 

.408268 

.409097 

.408268 

.408268 

.1 

.405724 

.406315 

.405754 

.405810 

.?}- 

.403534 

.403534 

.403533 

.403645 

.3 

.401648 

.402007 

.401628 

.401749 

.4 

.400026 

.400481 

.399979 

.400099 

.5<1 

.398629 

.398954 

.398629 

.398672 

.6L 

.397427 

.397427 

.397438 

.397445 

.7 

.396393 

.396590 

.396400 

.396395 

.8 

.395502 

.395752 

.395494 

.395499 

.9 

.394736 

.394914 

.394719 

.394733 

XLQC 

.394076 

.394076 

.394076 

.394076 

L  = 

knots  for  linear  spline,  Q  =  knots 

for  quadratic 

spl ine, 

C 

=  knots  for  cubic  s 

pi  i  ne 
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Fig.    1.     Plot   of  o+(.49-a)exp(-B(15t+17))   on   [-1,1]  with   a  =   .39  and 
0   =   .10  for   Example  2.2 


T 
G.475^v 

i   \ 

i 

Q.  425-1 


0. 400- 


-L.0 


-0.5 


0.0 


L.O 


Fig.  2.  Plot  of  Linear  Spline  Interpolate  to  a+( .49-a)exp (-S( 15t+17) ) 
with  Six  Evenly  Spaced  Knots  on  [-1,1]  for  Example  2.2 
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Q .  475-] 
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0.425- 


0.400- 
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-1.0 


-0.5 
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Fig.  3    Plot  of  Quadratic  Spline  Interpolate  to  a+( .49-a)exp(-B(15t+17) 
with  Five  Evenly  Spaced  Knots  on  [-1,1]  for  Example  2.2 


T 

0.  U7C1- 


0.450- 


0. 425- 


0. 400  J 


0.  375-^ 


1.0 


■0.5 


0 .  0 
X 


0.5 


Fig.  4.     Plot  of  Cubic  Spline  Interpolate  to  a+(  .49-a)exp(-3(15t+17) 
with  Three  Evenly  Spaced  Knots  on  [-1,1]  for  Example  2.2 


1,0 
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2.5  Discussion 


The  proposed  optimal  design  is  obtained  from  the  interpolating 
polynomials,  so  the  efficiency  of  the  design  depends  on  the  accuracy  of 
the  interpolations.  In  Example  2.2  we  compared  the  designs  based  on  the 
linear,  quadratic  and  cubic  splines  to  the  Box-Lucas  design.  Table  4 
summarizes  the  design  comparisons,  by  showing  that  the  optimal  designs 
based  on  linear  and  quadratic  splines  have  higher  efficiencies  than  the 
3ox-Lucas  design  but  the  optimal  design  based  on  cubic  spline  has  lower 
efficiency  than  the  Box -Lucas  design.  This  may  be  explained  by  the 
accuracies  of  the  interpolations;  that  is,  the  upper  bound  of  the  error 
of  interpolation  using  cubic  spline  is  larger  than  those  of  the  linear 
and  quadratic  splines  (see  Eq .  2.82).  We  also  note  here  that  the  design 
based  on  the  linear  spline  has  slightly  higher  efficiency  than  the 
design  based  on  the  quadratic  spline  when  comparing  both  of  them  to  the 
Box-Lucas  design.  Based  on  the  efficiencies  of  the  designs  our 
recommendation  for  Example  2.2  is  to  choose  the  design  that  is  based  on 
the  linear  spline.  The  number  of  design  points  (n)  for  the  linear, 
quadratic  and  cubic  splines  in  this  example  are  6,  6,  and  5, 
respectively.  This  number  seems  to  be  noni ncreasing  as  the  degree  of 
the  spline  interpolating  polynomial  increases.  This  is  not  always  true 
since  the  number  of  design  points  depends  on  m  and  h  (n  -  m+h+1)  where  m 
is  the  degree  of  spline  interpolating  polynomial  and  h  is  the  number  of 
interior  knots.  The  value  of  h  can  increase,  decrease  or  stay  the  same 
when  m  increases.  Hence,  there  is  the  possibility  of  increasing  the 
number  of  design  points  with  an  increase  in  the  degree  of  the  spline 
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interpolating  polynomial.  The  optimal  design  is  also  affected  by  the 
position  of  the  design  points  in  the  experimental  region.  A  problem 
arises  of  how  the  design  points  should  be  allocated,  especially  for  the 
design  that  is  based  on  spline  interpolating  polynomial.  Sox  and  Draper 
(1975)  point  out  that  "the  design  should  generate  a  satisfactory 
distribution  of  information  throughout  the  region  of  interest."  For 
this  reason,  for  the  optimal  design  based  on  spline  interpolating 
polynomials,  we  would  choose  our  design  so  that  there  is  at  least  one 
design  point  in  each  subinterval.  One  other  problem  is  that  when  the 
number  of  design  points  becomes  large  (n  >  10),  it  is  almost  impossible 
to  calculate  the  design  points  through  the  Controlled  Random  Search 
Procedure.  This  was  the  case  when  we  tried  to  calculate  the  design 
points  for  the  design  that  is  based  on  Lagrange  interpolating  polynomial 
with  two  input  variables  (Example  2.1).  It  is  hoped  that  a  sequential 
procedure,  in  which  points  are  added  one  at  a  time,  can  be  used  to  solve 
the  problem  of  optimization  when  the  number  of  design  points  is  large. 
This  is  subject  to  future  research. 


CHAPTER  THREE 

CONFIDENCE  REGIONS  FOR  THE  PARAMETERS 

IN  NONLINEAR  MODELS 


3.1  Introduction 


In  this  chapter  we  present  a  method  for  the  construction  of  a 
confidence  region  on  a  vector  of  parameters  in  a  nonlinear  model. 
Consider  the  same  model  as  in  (1.5) 

y  =  f(x;  8)  +  e   ,  (3.1) 

where  y  is  the  measured  or  observed  response,  f(x_;  _8_)  is  the  true  mean 
response,  _x_  =  (x^,  X2,  ...,  \)'  is  a  vector  of  k  input  variables,  Q_  = 
(8  8  ...,  8  )'  is  a  vector  of  p  unknown  parameters,  and  e  is  a 
random  error.  If  there  are  n  observations  (n  >  p)  of  the  form  yi,  y?, 
...,  yn  with  y.j  being  obtained  at  the  design  point  x.-j  =  (x-ji>  *i?>  •••, 
x-j^)',  i  =  1,  2,  ...,  n,  then  model  (3.1)  can  be  written  in  the  form 

y1  =  f(x.i  >   1)  +  £i   •  i  =  1,  2,  ....  n   .        (3.2) 
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We  assume  that  the  random  errors   e.  ,  e„,  ....  e   are  independent! v 

id  n         r      J 

p        p 

distributed  as  N(0,  a  ),  with  ac   unknown.   We  define  the  n  x  k  design 

matrix  as  D_  =  [x^],  i  =  1,  2,  ...,  n;  j  =  1,  2,  ...,  k .  The  i -th  row 
Xj  of  this  matrix  with  elements  x^,  x^p,  ...,  x^  provides  the  levels 
of  the  k  variables  at  which  the  response  is  to  be  observed  at  the  i  -th 
run.  Throughout  our  development,  we  shall  assume  that  replications  can 
be  obtained  at  each  design  point.  These  replications  will  be  used  to 
obtain  an  estimate  of  a2  (Draper  and  Smith  (1981)).  In  the  following 
sections,  we  consider  the  problem  of  constructing  a  confidence  region  on 
the  vector  of  parameters,  _8_. 

3.2  Method  for  Obtaining  Confidence  Region 
on  Parameters 

Consider  the  model  given  in  (3.1).  For  a  given  value  of  j<_,  f(x;  8) 
is  a  function  of  _8.  When  each  element  in  x_  takes  on  n  values  in  the 
experimental  region  R,  we  end  up  having  n  functions  of  6.  To  construct 
a  confidence  region  on  the  parameters  _6_,  we  shall  proceed  as  follows. 

Step  One.  For  n  >  p,  we  select  q  =  (n)  sets  of  points  each 
consisting  of  p  design  points.  Let  W  consist  of  the  collection  of  the  q 

sets  of  p  points.  Let  xj  ' )0  '     denote  the  p  points  in  the  1-th 

set  (1  =  1,  2,  ...,  q). 

Step  Two.   For  the  1-th  set  of  points  in  W,  1  =  1,  2,  ...,  q,  we 
obtain  a  set  of  simultaneous  confidence  intervals  on  f(x:  ';  8),  ..., 
f(/  '■   8_).   Let  a^1)  be  the  cartesian  product  of  these  intervals, 
that  is,  fn  '  is  a  rectangular  confidence  region  on  the  vector 


,'f(x 


:i 


..,  f(x 

-p 


(1 


Then  the  confidence  intervals  will  be 


-1  '  -' 
chosen  so  that  the  confidence  coefficient  for  this  region  is  equal  to 

1  -  a. 

Step  Three.  Let 

_(1 


'1 


f(x{1};  1) 


n^  =f(x^;  9)   , 


We  shall  solve  the  system  of  the  above  equations  for  9i,  . ..,  9  in 

terms  of  nj  ,  ...,  n^  '  (this  can  be  achieved  by  assuming  certain 

conditions  on  the  functions   ffx!1^;  9),  ...,  f(x^;  9),  (1  =  1,  2, 

— 1   —       *— ■ p   —   v       ' 


.  q 


Step  Four.   We  obtain  a  rectangular  confidence  region  on  a  vector 
with  a  confidence  coefficient  >  1  -  a.  This  region, 


1=   («!.  ••-,  % 

denoted  by  n1',  is  based  on  the  region  tt^>. 

Step  Five.    Using  Bonferroni's  inequality,  the  intersection, 

q   (1) 
O  rv     provides  a  confidence  region  on  the  vector  9  with  a 

1=1 
confidence  coefficient  >  1  -  qct.   This  intersection  can  be  used  to 

obtain  simultaneous  confidence  intervals  on  any  set  of  functions  of 

9i>  •••.  v 

We  shall  now  give  details  of  Steps  One  through  Five. 


3.3  Confidence  Region  on  _9 
In  this  section  we  first  construct  the  rectangular  confidence 
region  fiC)  on  each  of  a  set  of  p  functions  of  the  form  f  (xP  ' ;   9), 


69 


i=l,  2,  ...,  p,  where  xj  ',  ...,  x^  '  is  a  set  of  points  in  W, 
1  -  lj  2,  ...,  q.  The  p-variate  student's  t -distribution  defined  by 
Ounnett  and  Sobel  (1954)  will  be  used  to  obtain  these  confidence 
regions.  For  the  sake  of  simplicity,  in  this  section  we  shall  drop  the 
superscript  (1)  and  just  use  Xj . 

Consider  the  model  given  in  (3.2)  and  suppose 

yll'  VIZ'    ••*'  yini  are  nl  rePeat  observations  at  xj 
y21'  y22'  "*'  y2n9  are  n2  raPeat  observations  at  jo? 


ynl»  yn2»  •••>  ynn  are  np  repeat  observations  at  x^ 


where  yiu  is  the  u-th  observation  at  the  i -th  design  point.  Let  us 
assume  the  yi(J  are  independently  normal  distributed  with  mean  f (x_-  ;  9_) 
and  common  unknown  variance  a2  (i  =  1,  2,  ...,  n;  u  =  1,  2,  ...,  n-)- 
Draper  and  Smith  (1981)  pointed  out  that  to  estimate  a2  we  can  use  the 
pooled  variance  of  the  form 


"i      -  ? 

~       .1,      £  (y.  -y.) 
2  _  i=l  u  =  l  VJiu  J -[' 

1=1  ^i"1) 

(3  3) 
(nrl)s2+(n2-l)s^...+(nn-l)s2 
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where  s-   is  the  sample  variance  of  the  n1-(>  2)  observed  responses  at 

the  i-th  design  point,  y.  is  the  average  of  the  values  of  the  observed 

response  at  the  i-th  design  point,  and  n  is  the  number  of  distinct 
design  points.  Here  s'  has  m  =  •  Mn.-l)  degrees  of  freedom. 

Let 


*1 


y--f(x.; 


.i=l,2, 


,  p. 


3.41 


We  note  here  that  y\   is  normally  distributed  with  mean  f  (x, :  _6)  and 

o 

variance  ac/n^  .      If  we  write 


/n7  fy.-f (x. ;    6) 


i    =   1,    2, 


,   p, 


'3.5] 


then  z.j,  i  =  1,  2,  ...,  p,  has  a  p-variate  normal  distribution  with  mean 
0  and  common  unknown  variance  a2  and  Zi,  ,..,  z     are  independent  by  the 

assumption  regarding  the  random  errors  stated  in  Section  3.1.   Also, 

2  2 

ms  /a   has  a  chi  -square  distribution,  independent  of  the  z- ,  with  m 

degrees  of  freedom.  By  substituting  zi  into  (3.4)  we  obtain 


li   =T 


i  =  1,  2, 


3.6) 


The  joint  distribution  of  ti ,  i  =  1,  2,  ...,  p,  of  the  form  (3.6),  is  a 
central  p-variate  t-di stribution  with  m  degrees  of  freedom  and  has  the 
density   function 


71 


v   1  '     p' 


r(i(m+P; 


2V' 


m+p, 


-,  1   +  -  .  Z,    t. 

(mn)P/2rC^)         m  1=1     1 


,-°°<t.  <°°,   i=l,    . . . ,   p, 

(3.7) 


where  r(    )    is   the  Gamma   function    (see   Ounnett   and   Sobel    (1954)). 

Simultaneous  confidence  intervals  on  f  (x_.  ;  _9 ) ,  i  =  1,  2,  ....  p, 
with  confidence  coefficient  1-a  can  be  set  up  by  finding  p  positive  real 
numbers,   say   a^   a?,    -..,   a   ,   such  that 


P(|t1|<a1,|t2|<a2>. 


a       a     . 

_P       P-l 


K\<*J  -  /    / 


-a     -a     , 
P       P-l 


*   1   -  a 


J     h(tlf....t    )dtr..dt 


■3.8) 


3y    substituting    the   ti    given    in    (3.4),    formula    (3.8)    can   be   written   in 
the  equivalent   form 


p(yi-»l7=^(xi:  D<y1+a1— -,...,  y     a  -^<f(x    ;  e)<y  +a  -1-)  =  1-a 


3.9) 


which  determines  100(1  -  a)%  simultaneous  confidence  intervals  on 
f(x_1;  9.)  i  •••,  ?(**'  H)'  Let  n  denote  the  cartesian  product  of  these 
intervals.   The  region  fi  is  a  rectangular  confidence  region  on  the 
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vector  (f(x_i5  £) ,  •••,  f(_*n;  jO ) '  wl"tn  a  confidence  coefficient  1  -  a. 

Let  n.j  =  f  (x_.  ;  6_) ,  i  =  1,  2,  ...,  p.  ni  i  s  a  function  of  x_j  and 
9.  For  a  given  value  of  x* ,  f  (x_-  ;  8_)  is  a  function  of  _9  and  we  write 
f  (x_.  ;  8_)  as  f-j  (_6_) .  Thus  we  can  express  ni ,  i  =  1,  2,  . .  • ,  p,  as 


nL  =f1(61,  82. 
n2  =  f2^9l'  92' 


V  ■  fid) 


n  «f  (8,,  8 ,  8  )  «fj 

p    pv  1'  2'      p'    p 


Using  the  expressions  in  (3.10)  we  can  write  (3.9)  as 


3.10) 


P((nr  n2, 


...,  n)efl)al-a 


3.11) 


If  the  Jacobian 


3(fp 


'.  fJ 


3(0 


£y  *  0,   then  equations  (3.10)  can  be 


r 


uniquely  solved   for  8,,    ...,   8Q  in  terms   of  m,    •..,   n»  and  we  obtain 


*2  =  92^ n 1 »   n2! 


.   \ 


>    *V 


'3.12] 


»p  =  gp(n1>  n2,    ...,   np) 


We  let  L.  =  min  g.  (n,,  ...,  n  )  and  U.  =  max  g.  (n, 
i    Q  yi   1       p        i    fl   yi   1 

from  (3.11)  and  (3.12)  we  can  conclude  that 


,  n  ),  then 
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P(Ll  <  8X  <  U^  ....  Ip  <  9p  <Up)  >  1  -  a 


3.13) 


which  gives  a  rectangular  confidence  region  on  _9_  with  a  confidence 
coefficient  >  1  -  a.  We  denote  this  region  by  r.  Equation  (3.13)  can 
be  written  in  the  alternative  form 


P((81,  02,  ...,  8  )  e  T)  >  1  -  o 


:3.h) 


The  above  method  can  be  applied  to  every  set  in  W  consisting  of  p 
points.  We  shall  denote  the  confidence  region  obtained  as  in  (3.14)  for 


the  various  q  sets  in  8  by  r 


i: 


We  thus  have  the  following 


q  confidence  regions  on  8_  each  having  a  confidence  coefficient  >  1  -  a, 


p((8lf  e2,  ...,  8p)  £  r(1))  >  l  -  a 
P((e1,  e2,  , 


9   e  r 

P1 


>   1 


'3.15) 


P((8r  e2,  ...,  8  )  E  r(^})  >  1  - 


We  now  combine  the  information  from  these  q  confidence  regions  using 
Bonferroni  's  inequality  (see,  for  example,  Meter  and  Wasserman  (1974,  p. 
146)  or  Bickel  and  Doksum  (1977,  p.  439)).  According  to  this  inequality 
and  from  (3.15),  we  have 


P((er 


2' 


q 

=  n 

1=1 


r(l) 


>  1  -  qa 


:3.iR) 
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which  gives  a  confidence  region  on  _9  with  a  confidence  coefficient 
>  1  -  qo.  We  can  write  (3.16)  in  the  equivalent  form 

P(L*  <  8X  <  U*.  ...,  L*  <  6p  <  U*)  >  1  -qa   ,      (3.17) 

L*  =  maxdj1^  ...,  LJq)),  U*  =  min(uj1),  ....  u{^),  1  -  1,  2,  ...,  p, 
where  l{1},  uj]),  i  -  1,  2,  ...,  p,  1=1,  2,  ...,  q,  are  the  lower 
and  upper  bounds,  respectively,  on  8i  obtained  as  in  (3.13)  for  the  1-th 
set  of  points  in  W.  We  note  from  formula  (3.17)  that  the  intervals 
[L*,  U*],  i=l,  ...,  p,  form  a  set  of  simultaneous  confidence  intervals 
on  9i-  ■••>  9p,  respectively,  with  a  joint  confidence  coefficient 
greater  than  or  equal  to  1  -  qa. 


3 . 4  Examples 

We  now  present  two  examples  to  illustrate  the  technique  for 
obtaining  the  confidence  regions  of  parameters  in  nonlinear  models. 

Example  3.1.  Consider  the  two  parameter  exponential  growth  model 
used  by  Gallant  (1975b)  of  the  form 


ex 

y  =  V   +  e  •  (3.i8) 


A  total  of  n  =  11  design  points  are  chosen  from  R  =  {x:  0  <  x  <  l}. 
These  points  are  0,  .1,  .2,  ...,  1.   Each  point  is  replicated  five 
times.   The  parameter  space  is  taken  as  0  =  {(8,,  8„):  6  >  0,  8o  >  0}. 
Observed  response  values  are  generated  using  the  parameter  values, 
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8^  =  .62,  &2    =  *5,  and  a^  =  .04.   These  values  are  shown  in  Table  5. 

9 

The  mean  (y.  )  and  the  variance  (s.)  of  the  observed  responses  at  the 

i-th  design  point  for  i=l,  ...,  11,  are   listed  in  Table  7.  The  pooled 

variance  s  ,  can  be  calculated  by  using  the  formula  in  (3.3)  and  is 

equal  to  .0332818  with  m  =  44  degrees  of  freedom. 

Since  n  =  11   and  p,  the  number  of  parameters  is  equal  to  2,  then 

q  =  (  _)   =55.   We  thus  have  55  pairs  of  the  functions  of  the  form 

f(xP);_e),   1  ■  1,  2;  1  -  1,  2,  ...,  55.   For  each  pair,  the 

simultaneous  confidence  interval  on  f(xj   ;  _9 ) ,  f(xi  ,  _9 ) ,  can     be 

obtained  as  follows.   First  we  use  the  bi-variate  t-distribution  given 

in  (3.7)  with  p  =  2  to  calculate  the  probability  described  in  (3.8).  We 

assume  that  the  limits  of  integration  a^,  a2  are  equal  to  a,  say,  so 

that  Equation  (3.8)  can  be  expressed  as 

a  a 
P(  |t1|  <  a,  |t2|  <  a)  =  /  /  h(tr  t2)  dt,dt2  »  1  -  o 

-a  -a 

(3.19) 


where  h(t1,  t2)  can  be  obtained  from  (3.7)  with  p  =  2  and  m  =  44.  To 
calculate  the  probability  value  in  (3.19)  we  fix  a  and  calculate  the 
integrals.  This  can  be  carried  out  by  using  a  computer  subprogram 
called  DMLIN  from  IMSL  (1982).  Suppose  we  would  like  the  confidence 
coefficient  1  -  qa  in  (3.16)  to  be  equal  to  .94.  We  set 
1  -  (  2)a  =  .94  which  implies  that  a  :  .001.  3y  computing  the 
integral  in  (3.19)  for  several  assignments  of  a  we  obtain  a  =  3.75  which 
corresponds  to  a  =  .001057813  and  1  -  (  ha  =    .941319999.   To  find 
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Table  6 
Observed  Responses  and  Design  Points  for  Example  3.1 


Observed 
Responses 

Design 
Points 

Observed 
Responses 

Design 
Points 

Observed 
Responses 

Design 
Points 

(y) 

(x) 

(y) 

(x) 

(y) 

(x) 

.84025 

0 

.65352 

.3 

.82359 

.7 

.84356 

0 

.85087 

.4 

.63941 

.7 

.66729 

0 

.63028 

.4 

.63028 

.7 

.75769 

0 

.56506 

.4 

.94460 

.8 

.67426 

0 

.81007 

.4 

.71885 

.8 

.53177 

.1 

.89356 

.4 

.89828 

.8 

.62303 

.1 

1.10607 

.5 

1.00649 

.8 

.75965 

.1 

1.09510 

.5 

.74048 

.8 

.49493 

.1 

.64805 

.5 

.96815 

.9 

.70555 

.1 

.88696 

.5 

1.08855 

.9 

.77509 

.2 

.85727 

.5 

.66642 

.9 

.61030 

.2 

1.31870 

.6 

.79255 

.9 

.79161 

.2 

.73195 

.6 

.88114 

.9 

.58948 

.2 

.85361 

.6 

1.0697 

1 

.74782 

.2 

.52289 

.6 

.56959 

1 

.46908 

.3 

1.05243 

.6 

.91519 

1 

.82377 

.3 

1.02124 

.7 

1.31447 

1 

.91124 

.3 

1.05300 

.7 

.84882 

1 

.66913 

.3 
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Table  7 
Means  and  Variances  of  Observed  Responses 
at  Distinct  Design  Points  for  Example  3.1 


Design  Points 

(Xi) 


Observed  Responses 


Mean  (y. 


Variance  fs' 


0 
.1 
.2 
.3 
.4 
.5 
.6 
.7 
.8 
.9 
1 


.75661 
.62299 
.70286 

.70635 
.74997 
.91869 
.89592 
.83350 
.86174 
.87936 
.94355 


.00733 
.01255 
.00913 
.02893 
.02073 
.03606 
.09277 
.04062 
.01607 
.02613 
.07578 
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the       100(1   -  a)%  *         99.9%       simultaneous       confidence       intervals       on 
f{x[]>;   9_),      and      f(x^;_6),   1    =    1,    2,    ...,    55,    we    consider    Equation 
(3.9)    with    a]_    =   a2   =   a   =   3.75,    rii    =   r\2   =    5.      These   intervals    are   given 
in   Table  8. 

We  now  find  the  confidence  region  on  (8,,  62).  For  the  sake  of 
simplicity,  in  this  example  we  shall  drop  the  superscript  (1)  and  just 
use  x^  and  Xg-     We  assume  that  Xi    <  Xo  and   let 


9  x 
nx   =  8ie  f(xx;    0) 


!3.20) 


V2 
n2  ,  9ie  f(x   :   9) 


Suppose  Cj,  di  are  the  lower  and  upper  limits  of  the  confidence  interval 
on  f(x^;  9_)  and  c^,  d2  are  the  lower  and  upper  limits  of  the  confidence 
interval  on  f(x2;  9_) .  We  solve  the  two  equations  in  (3.20)  for  9.,  8 


and  obtain 


A  -*        A  s\  A  r\  "~A   1 

n2  nx 


In  n^-ln  n. 


3.21! 


%2~*1 


where  In  is  the  natural  logarithm.  From  (3.21)  we  obtain  by  minimizing 
and  maximizing  each  of  8,  and  92  with  respect  to  n  over  the  rectangular 
region  [cj_,  dj]  x  [c2,  d2] 
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1  _^2__ 

A   -I    ~A  ^  X  q  — X   -1 

lmin       a2  cl 


Xl  X2 


^1    "An  Xrt"~X^ 

'lmax  =  c2  dl 


2mi  n 


In  c„-ln   d. 


2max 


In  dp-ln   c, 
x2"xl 


Values  of  e^,  9lmax  and  e2min,  92max  for  each  of  the  55  pairs  of 
simultaneous  confidence  intervals  on  ni ,  n2  are  given  in  Table  8.  From 
Table  8,  if  we  take  the  largest  of  the  lower  bounds  on  8j  and  on  92  as 
well  as  the  smallest  of  the  upper  bounds  on  91  and  on  9o,  we  obtain  the 
following  rectangular  confidence  region  on  Q_  which  has  a  confidence 
coefficient  of  at  least  94.18  percent: 

8j  e  [.45065,  .96861] 

62  e  [-.41816,   1.01981] 


For  comparison  purposes  we  use  the  data  in  Table  6  to  construct 
simultaneous  confidence  intervals  on  8j  and  e2  by  using  each  of  three 
other  known  methods.  These  methods  are: 
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Table  8    m       m 
99.9  Percent  Confidence  Regions  on  [f(x1lu;  _9 ) ,  f(x2       >   ±)J 

1  =1,2,  ,..,  55  and  at  least 

99.9  Percent  Confidence  Regions  on  [8i,  82]  for  Exanple  3.1 

Design  Points  Confidence  Region  of  Confidence  Region  of 

x|])      XJ)  [fOcJ13;  e),  fOcCl).  e)]  fy,  e2] 


0  .1  [.45065, 1.06257], [.31703, .92894]  [.45065,1.06257],[-12.094,7.2335] 

0  .2  [.45065, 1.06257], [.3969,1. 00882]  [.45065,1.06257],  [-4.9238,4.02918] 

0  .3  [.45065, 1.06257] ,[.  40039, 1.0123]  [.45065, 1.06257],  [-3.2533,2.69762] 

0  .4  [.45065, 1.06257], [.44401 ,1.05592]  [.45065, 1.06257],[-2. 1815, 2. 12369] 

0  .5  [.45065,1. 06257]  ,[.61273, 1.22465]  [.45065, 1.06257],  [-1.101, 1.99942] 

0  .6  [.45065, 1.06257], [.58996, 1.02187]  [.45065, 1.06257], [-.98065, 1.6349] 

0  .7  [.45065,1. 06257], [.  52755, 1.13946]  [.45065, 1.06257],  [-1.0003,1.32516] 

0  .8  [.45065, 1.06257], [.55578,1. 1677]  [.45065,1.06257], [-.81008,1. 19011] 

0  .9  [.45065, 1.06257], [.57341, 1.18532]  [.45065, 1.06257], [-.68539, 1.07452] 

0  1  [.45065,1.06257], [.63759, 1.24951]  [.45065, 1.06257], [-.51074, 1.01981] 

.2  [.31703,.  92894], [.3969, 1.00882]  [.09963,2. 17  417], [-8.5036, 11.5754] 

.3  [.31703,. 92894],[.40O39, 1.0123]  [.17742, 1.41495],  [-4.208, 5.80496] 

.4  [.31703, .92894], [.44401, 1.05592]  [.21229, 1.1881], [-2. 4607, 4.0106] 

.5  [.31703,. 92394], [.61273, 1.22465]  [.22614,1.03079],  [-1.0403, 3.37354] 

.6  [.31703,  .928^], [.58996, 1.02137]  [.24236, 1.01724],  [-.90799,2.66529] 

.7  [.31703,. 92394],[.52755, 1.13946]  [.25615, 1.02081], [-.94302, 2. 1322] 

.8  [.31703,. 92894], [.55578,1. 1677]  [.26315,. 99967], [-.73381, 1.86256] 
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Table  8— continued. 

.1  .9  [.31703,. 92894], [.  57341,1.18532]  [.26885,.  98669]  ,[-.60307, 1.64347] 

**  * 

.1  1  [.31703,. 92894], [.63759, 1.24951]  [.27222,. 96861], [-.41816,1.5239] 

.2  .3  [  .3969, 1 .00882]  ,C  .40039, 1 .0123]  [.06101,6. 40426], [-9. 2409,9.36293] 

.2  .4  [.3969,1.00882], [.44401, 1.05592]  [.14919, 2.29209],[-4. 1034,4.8924] 

.2  .5  [.3969, 1.00882], [.61273, 1.22465]  [. 18727, 1.40661] ,[-1 .662,3.75572] 

.2  .6  [.3969,1.00882],  [.58996,1.02187]  [.22809, 1.31919], [-1.3412, 2.76986] 

,2  .7  [.3969, 1.00882], [.52755, 1.13946]  [.  26030, 1.30747], [-1.2966, 2.10924] 

.2  .8  [.3969,1.00882], [.55578,1. 1577]  [.27699, 1.23059], [-.99359,1.79849] 

.2  .9  [.3969, 1.00882], [.57341, 1.18532]  [.29035, 1.18553], [-.80706, 1.56296] 

.2  1  [.3969,1.00832], [.63759,1.24951]  [.29797,1.13143], [-.57353, 1.43352] 

.3  .4  [.40039, 1.0123], [.44401, 1.05592]  [.02133, 11. 9968], [-8.2414,9. 69731] 

.3  .5  [.40039, 1.0123],  [.61273, 1.2246  5]  [.07435, 2. 14967],  [-2.5103, 5.58933] 

.3  .6  [.40039, 1.0123], [.58996, 1.02137]  [. 13339, 1.737],[-1. 79998,3.66398] 

.3  .7  [.40039, 1.0123],  [.52755, 1.13916]  [.18274, 1.65044],  [-1.6294, 2. 61467] 

.3  .8  [.40039, 1.0123] ,[.  55578,1. 1677]  [.21066,1. 45062], [-1.1992, 2. 14069] 

.3  .9  [.40003, 1.0123], [.57341, 1.13532]  [.23271 ,1.34504], [-.94732, 1.80883] 

.3  1  [.  40039, 1.0123],[.  63759, 1.24951]  [.24535, 1.23411], [-.6604, 1.62581] 

.4  .5  [.44401 ,1.05592], [.61273, 1.2246 5]  [.00767,9.31279], [-5.4424, 10. 1456] 

.4  .6  [.44401, 1.05592], [.58996, 1.02187]  [.06059,3.38264], [-2. 9106,4.97893] 

.4  .7  [.44401, 1.05592], [.52755, 1.13946]  [.12637,2.66356], [-2.3131, 3.14154] 

.4  .8  [.44401, 1.05592], [.55578,1. 1677]  [. 16883, 2.00614], [-1.6045,2.41735] 

.4  .9  [.44401, 1.05592], [.57341, 1.18532]  [.20241, 1.72(396],  [-1.2212 ,1.96384] 

.4  1  [.44401, 1.05592],  [.63759, 1.24951]  [.  22276, 1 .47806] ,[-  ..84378,1. 72443] 

.5  .6  [.61273, 1.22465], [.58996, 1.02187]  [.02110,47.2015], [-7.3035, 6. 73707] 
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Table  8— continued. 

.5  .7  [.61273, 1.22465], [.52755,1. 13941]  [.12993, 10.0551], [-4.21(18,3. 10191] 

.5  .8  [.61273, 1.22465], [.55578, 1.1677]  [.20917,4. 56934], [-2.6334, 2. 14953] 

.5  .9  [.61273, 1.22465], [.57341, 1.18532]  [.26858, 3. 1619], [-1.897, 1.64959] 

.5  1  [.61273, 1.22465], [.63759, 1.24951]  [.30047, 2. 35221], [-1.3054,1.42516] 

.6  .7  [.58996,1.02187],[. 52755,1. 13946]  [.01136, 168.052], [-8.234,6. 58257] 

.6  .8  [.58996, 1.02187], [.55578, 1.1677]  [.07608,12. 154], [-3.8563, 3.41368] 

.6  .9  [.58996,1.02198],[. 57341, 1.18532]  [.14615,5.28023], [-2.4663,2.32571] 

.6  1  [.58996, 1.02187], [.53759,1.24951]  [.19140, 3. 11048],  [-1. 5848, 1.87614] 

.7  .8  [.  52755, 1.13946] ,[.  55578,1. 1677]  [.00203,173.484],[-7. 1793, 7.94551] 

.7  .9  [.52755,1.13946],[.57341, 1.18532]  [.03103,12.6047],[-3.4336,4.04765] 

.7  1  [.52755,1. 13946], [.63759,1.24951]  [.07055,4.41628], [-1.9353,2.87423] 

.8  .9  [.55578,1. 1677],[.57341, 1.18532]  [.00129,345.366],[-7.112,7.57389] 

.8  1  [.55578, 1.1677],  [.63759, 1.24951]  [.02176, 13. 1361], [-3. 0254,4.05065] 

.9  1  [.57341,1. 18532],  [.63759, 1.24951]  [.00052, 314.366],[-6.2O06,7.78915] 


*  is   the   largest   value  of  the  lower  bound. 
**  is   the  smallest   value  of   the  upper  bound. 
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1)  The  likelihood  ratio  method 

2)  The  method  based  on  the  asymptotic  normality  of  the  least 
squares  estimator 

3)  Hartley 's  method 
(see  Section  1 .3) . 

The  boundaries  of  the  confidence  regions  on  9i  and  8?  obtained  from 
the  above  three  methods  are  ellipses.  In  order  to  be  able  to  compare 
our  simultaneous  confidence  intervals  on  6,  and  82,  which  we  have  found 
earlier,  to  those  obtained  by  each  of  the  above  methods,  we  project  the 
ellipses  onto  the  8,  and  82  axes.  The  resulting  intervals  form 
simultaneous,  but  conservative,  confidence  intervals  on  8,  and  8?.  The 
results  are  given  in  Table  9. 


Table  9 
Simultaneous  Confidence  Intervals  on  8.  and  8- 

with  a  Confidence  Coefficient  >  .9418  Obtained  from  Projecting  the 

94.13  percent  Confidence  Regions  on  (8.,  6?)  for  Example  3.1 


Method 


Lower 
Limit 


Upper 
Limit 


Lower 
Limit 


Likelihood  Ratio  .585    .790    .09 

Asymptotic  Normality  of  the  LS  Estimator   .602    .768    .129 
Hartley  .555    .820    .02 


Upper 
Limit 

.55 
.506 

.63 
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We  note  that  the  proposed  method  provides  wider  confidence 
intervals  than  the  other  three  methods.  This  is  attributed  to  the 
conservative  nature  of  the  proposed  method  caused  by  the  use  of  the 
confidence  intervals  in  (3.13),  which  are  already  conservative,  followed 
by  the  use  of  Bonferroni's  inequality  which  leads  yet  to  another  set  of 
conservative  intervals  as  in  (3.17). 

Example  3.2.  In  this  second  example  we  construct  simultaneous 
confidence  intervals  on  9.,  9„  for  a  nonlinear  model  with  two  input 
variables.  Consider  the  model  used  by  Draper  and  Smith  (1981)  for  a 
certain  chemical  reaction  which  can  be  described  by  the  nonlinear  model 

y  =  expf-9^  exp[-62(l/x2  -  1/620)]}  +  e   ,      (3. 22) 


where  y  is  the  fraction  of  original  material  remaining,  x^_  is  the 
reaction  time  in  minutes,  and  x2  is  the  temperature  in  degrees  Kelvin. 
The  data  are  taken  from  page  520  of  Draper  and  Smith  (1981)  and  are 
given  in  Table  10.  The  data  consist  of  two  replications  at  x^  =  120, 
x2  =  612,  eleven  replications  at  x^  =  60,  x2  =  620,  and  three  replica- 

o 

tions  at  x1  =  30,  x2  =  639.  The  mean  (y.  )  and  variance  (s  )  of  the 
observed  responses  at  each  setting  of  x^  and  x2  are  given  in  Table  11. 

Before  we  proceed  any  further  we  would  like  to  test  to  make  sure 
that  the  variance  of  the  observed  responses  from  each  set  of  replica- 
tions are  not  significantly  different  because  in  our  method  we  use  the 
pooled  variance.  To  test  the  null  hypotheses  of  equal  variances, 
Bartlett's  (1937)  chi-square  test  is  considered  (see  also  Brownlee 
(1965,  p.  293)).  The  test  statistic  is 
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M  =  [(N  -  n)  In  s2  -  Z     (n.  -  1)  In  s2]/C   ,      (3.23) 
where 

and  s2  =   .00006222  is   the  pooled  variance  defined   in   (3.3), 

n  2   2   2 

N  =  Z     n.  =  16,   s,,  s9,  s^  are  in  Table  11  and  n  is  the  number  of 
i=1  i        i  i.       j 

distinct  design  points.  Under  the  null  hypothesis  M  has  an  approximate 
chi -square  distribution  with  n  -  1  =  2  degrees  of  freedom.  The  value  of 
M  in  (3.23)  is  equal  to  .988932.  The  value  of  chi -square  with  2  degrees 
of  freedom  and  the  significance  level  of  .05  is  equal  to  5.99  which 
exceeds  the  value  of  M.  Therefore,  we  do  not  reject  the  null  hypothesis 
of  equal  variances. 

We  have  n  =  3  distinct  sets  of  design  points  with  two  parameters. 
There  are  (2)  =  3  pairs  of  functions  of  the  form  f  (x.  V  ,x  )V ;  9_) , 
i  -  1,  2;  1  =  1,  2,  3.  To  obtain  a  confidence  region  on  _9  with  a 
confidence  coefficient  of  at  least  95  percent  we  need  a  to  satisfy 
1  -  L)a  =  .95,  thus  a  Z  .015.  3y  computing  the  integrals  in  (3.19) 
for  several  assignments  of  a  we  obtain  a  =  3.15  which  makes  1  -  o  = 
.984949706  and  1  -  [*)*  ■  .954849117.  The  98.49  percent  simultaneous 
confidence  interval  on  f  (x[J  ,x{^  ;  _6) ,  f  (xPJ  ,x^J  ;  _9) ,  1=1,2,  3, 
can  be  obtained  from  Equation  (3.9)  with  p  =  2,  m  =  13,  n,  =  2,  n?  =  11, 
n3  =  3  and  y . ,  s^  are  as  in  Table  11.  The  results  are  given  in 
Table  12. 

Again,  in  this  example  we  shall  drop  the  superscript  (1)  and  just 
use  Xj^,  x^»  *21  an<^  x2?" 
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Table  10 
Observed  Responses  and  Design  Points  for  Example  3.2 


Observed  Responses 

(y) 


.785 
.791 
.787 
.782 
.795 
.800 
.790 
.802 
.802 
.804 
.794 
.804 
.799 
.655 
.638 
.657 


Design  Points 

xl 

x2 

120 

612 

120 

612 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

60 

620 

30 

639 

30 

639 

30 

639 
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Table  11 

Means  and  Variances  of  Observed  Responses 

at  Each  Replication  for  Example  3.2 


Design  Points  Observed  Responses 

_  o 

X]_  x2  Mean    (y. )  Variance   (s." 


120 

612 

.78800 

.000018 

60 

620 

.79627 

.000054 

30 

639 

.65067 

.000124 

To   find  simultaneous  confidence  intervals   on   8.,   8_,let 


nl  =  f(xH'xi2;  D    =  exp{-91xnexp[-e2(l/x12  -   1/520)]; 
n2  =  f(x21,x22;  £)    =  exp{-6lX21exp[-e2(l/x22  -   1/620)]; 


3.24) 


Let  Cj,  d^  be  the  lower  and  upper  limits  of  the  confidence  interval  on 
f(x,,,x,-;  _e)  and  C2,d2  be  the  lower  and  upper  limits  of  the  confidence 
interval  on  fUoi'*??'  6_) .  We  solve  the  two  equations  in  (3.24)  for  9, 
and  8~  and  obtain 


Bj  =    (-In   nx)    exp[92(l/x12  -    l/620)]/xn 

(3.25) 

ln(-ln   ti   )    -   1  n{  -1  n   nJ    -   In   x..    +  In   x?, 
2  (l/x22-l/620)-(l/x12-l/620) 

We  now  find  the  maximum  and  minimum  values  of  8,  and  82  of  tne  form 
(3.25)  over  the  region  [c^d^  x  [c2,d2],  where  n1  e  [c^.d^]  and  n2  e 
[c2,d2].  The  results  are  also  in  Table  12.  The  largest  of  the  lower 
bounds  on  8j  and  92  and  the  smallest  of  the  upper  bounds  on  6,  and  9o 
from  Table  12  form  a  confidence  region  on  _9  with  a  confidence 
coefficient  of  at  least  95.48  percent  and  is  given  by 

81  e  [.00364,  .00392] 

92  e  [26558.7,  29617.4] 

For  comparison  purposes  we  use  the  data  in  Table  10  to  construct 
simultaneous  confidence  intervals  on  9  9  by  using  each  of  the  same 
three  methods  as  in  Example  3.1.  The  boundaries  of  these  confidence 
regions  are  ellipses.  The  simultaneous  confidence  intervals  on 
81  and  92  can  be  obtained  by  projecting  the  ellipses  onto  the 
61  and  92  axes.  These  confidence  intervals  are  conservative  and  are 
given  in  Table  13. 

Our  method  seems  to  do  well  when  comparing  the  ranges  of  the 
parameters  with  the  other  three  methods. 
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Table  12 
93.49  Percent  Confidence  Region  on 

Cf(xll  'x12  ;  -}'  f(x21  ,X22  ;  -)]>  ]    =  l'   Z>   3>  and  at  1east 
98.49  Percent  Confidence  Regions  on  [8^,  62]  for  Example  3.2 


Design  Points 

x(1))  (x(1)    x(1 
*12  ;  lx21  '    22 


Confidence  Region  on 


Cf(x{,1,^2!);9).f(xV,;,x-V;9)] 


<(1)x(1 

(21  'x22  '  -' 


Confidence  Region  on 
[er  e2] 


(120,612)  (60,620)    [.77043, .80557], [.78878, .80376]  [.00364, .00395], [24470.3,37284.9] 

! 120,  612)  (30,639)    [.77043,. 80557], [.63632,. 66501]  [.00334,. 00392], [26558. 7, 30762] 

;  60,  620)  (30,  639)     [.78878, .80376], [.63632, .66501]  [.00364,.00395],[25753.9,29617.4] 

*     is  the  largest  value  of  the  lower  baind 
**    is  the  smallest  value  of  the  upper  bound 


Table  13 

Simultaneous  Confidence  Intervals  on  8  and  9?  with  a  Confidence 

Coefficient  >  .9548  Obtained  from  Projecting  the 

95.48  Percent  Confidence  Regions  on  (8.,  8?)  for  Example  3.2 


Method 


'1 
Lower 

Limi  t 


Upper 
Limit 


Lower 
Limit 


Upper 
Limit 


Likel ihood  Ratio 

Asymptotic  Normality  of  the 
LS  Estimator 


.00361    .00389     26763      29393 
.00363    .00336     27006.3    29153.1 


Hartley 


.00356 


.00338 


27100 


30950 
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3.5  Discussion 


Our  confidence  region  on  the  parameters  is  conservative;  that  is, 
it  has  a  confidence  level  of  at  least  109(1  -  qa)  percent,  where 
q  =  [  ),  n  is  the  number  of  distinct  design  points  and  p  is  the  number 
of  parameters.  This  might  explain  the  reason  in  Example  3.1  that  our 
method  provides  wider  confidence  intervals  than  the  methods  of 
likelihood  ratio,  asymptotic  normality  of  the  least  squares  estimator 
and  Hartley.  The  advantage  of  the  proposed  method  is  that  it  does  not 
require  specifying  the  initial  estimates  of  the  parameters,  unlike  the 
likelihood  ratio  and  asymptotic  normality  methods.  Furthermore,  in 
comparison  with  Hartley's  method,  the  proposed  method  uses  the  true 
response  function  rather  than  the  linearized  one.  In  Example  3.2  our 
method  did  quite  well  by  comparison  with  the  other  three  methods. 

The  conservative  nature  of  the  proposed  method  means  that  it 
guarantees  a  percent  coverage  of  the  parameter  vector  which  is  no  less 
than  a  preset  value.  It  would,  therefore,  be  desirable  to  determine  how 
much  larger  is  the  true  percent  coverage  than  the  preset  value.  Such 
determination  can  help  reduce  the  size  of  the  confidence  region  obtained 
by  this  method.  Further  investigation  of  this  will  be  needed. 


CHAPTER  FOUR 
CONCLUSIONS  AND  RECOMMENDATIONS 


A  design  criterion  for  parameter  estimation  which  is  not  dependent 
on  the  parameter  values  in  a  nonlinear  model  as  well  as  a  method  for 
obtaining  a  confidence  region  on  the  parameter  vector  in  a  nonlinear 
model  have  been  proposed. 

In  Chapter  Two  approximations  to  the  true  mean  response  functions 
with  a  Lagrange  interpolating  polynomial  (for  the  case  of  two  input 
variables)  and  spline  interpolating  polynomials  (for  the  case  of  a 
single  input  variable)  were  constructed.  Each  approximation  was  used  to 
obtain  a  parameter-free  design  on  the  basis  of  a  modification  of  the 
D-optimality  criterion  which  was  proposed  by  Khun  (1982).  This  design 
is  found  to  depend  on  the  size  of  the  error  associated  with  the 
approximation  of  the  true  mean  response.  If  the  error  of  an 
approximation  is  large,  the  design  determined  on  the  basis  of  our 
procedure  may  be  inefficient. 

It  is  unfortunate  that,  at  the  present  stage,  no  computer  algorithm 
is  available  for  the  determination  of  the  optimal  design  based  on  the 
Lagrange  interpolating  polynomial  for  the  case  of  two  (or  more)  input 
variables  when  the  number  of  design  points  is  large.  It  is  hoped  that  a 
sequential  procedure  for  the  generation  of  design  points  can  reduce  the 
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computational   difficulty  to  a  manageable  level.    This  will   be 
investigated  in  future  research  (see  Example  2.1). 

The  reason  that  we  chose  spline  interpolating  polynomials  to 
approximate  the  mean  response  function  is  because  of  the  nice  properties 
of  the  splines  as  approximating  functions  (Wold  (1974)).  These 
functions  seem  to  be  extremely  suitable  for  interpolating  or 
approximating  functions;  to  quote  deBoor  and  Rice  (1968,  p.  1):  "Spline 
functions  are  the  most  successful  approximating  functions  in  use 
today.  They  combine  ease  of  handling  in  a  computer  with  great 
flexibility,  and  are  therefore  particularly  suited  for  the  approximation 
of  experimental  data." 

Our  attention  has  centered  on  the  construction  of  an  optimal  design 
based  on  spline  interpolating  polynomials  for  the  case  of  a  single  input 
variable.  Extensions  of  the  procedure  for  constructing  optimal  designs 
using  spline  interpolating  polynomials  to  more  than  one  input  variable 
cause  no  difficulties  in  terms  of  methodology.  However,  in  order  to 
calculate  the  design  point  locations,  we  may  have  the  same  difficulty  as 
in  Example  2.1  if  the  number  of  design  points  is  large. 

Our  design  criterion  was  developed  under  the  assumption  that  the 
values  of  all  the  parameters  are  of  equal  interest  to  the 
experimenter.  It  is  sometimes  the  case,  however,  that  the  experimenter 
is  more  interested  in  the  values  of  certain  parameters  than  the  values 
of  others.  The  extension  of  our  design  criterion  to  obtain  the  design 
of  experiments  for  subsets  of  parameters  requires  further  research. 
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In  Chapter  Three  a  method  for  the  construction  of  a  confidence 
region  for  the  parameters  in  nonlinear  model  was  developed.  Using  the 
confidence  intervals  on  the  mean  response  functions  and  the  Bonferroni's 
inequality,  a  conservative  confidence  region  for  the  parameters  was 
constructed.  Unlike  most  procedures,  this  procedure  does  not  require 
specifying  the  initial  estimates  of  the  parameters. 

As  we  mentioned  above,  our  method  yields  a  conservative  confidence 
region.  Further  research  needs  to  be  made  to  reduce  the  size  of  the 
confidence  region  without  affecting  its  confidence  coefficient.  This 
will  depend  on  gaining  some  knowledge  about  the  true  probability  of 
coverage  of  the  conservative  confidence  regions.  Furthermore,  an 
extension  of  our  method  to  construct  confidence  regions  on  subsets  of 
parameters  will  be  attempted. 


APPENDIX  1 
PROOF  OF  INEQUALITY  (2.78) 


Consider  model  (2.77).   The  (mj+D-th  partial  derivative  of  n(ti, 
t2;  _8_)  with  respect  to  ti  is 


m,       m.+l   m.+l 
^1+1)(t  ,  t  ;  8)  -  ('1}  V1*'9!   93bI   C2+92b2(t2+l)] 

1  "      L  [2+e1b1(t1+l)+92b2(t2+l)]mi+2 


and   the    (m2+l)-th   partial    derivative  of  n(tx,   t2;  _9)   with    respect  to  t? 


is 


^?+l  m,+l  m0+l 

(m2+1),  (-1)         (m-+l)!e.8,b1(t1+l)8/     b/ 

%  (t,,   t0:   8)   = ^ U_±_i £ £__ 

fc2         v    1       2     -;  m  +2 

C2+e1b1(t1+i)+e2b2(t2+i)3  < 


[A1.2! 


Upon  substituting  (Al.l)  and  (A1.2)  into  (2.51)  we  obtain 
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■Kt^l 
•Kt2<l 


!e(t1,  t2; 


2  -Kt,<l 


■1)         eie3bi(ti+1)e2       b2 

"~  m2+2    I 


"l        [2+e  b  (t1+l)+e  b?(t?+i)] 
-Kt2<l  :    L     *  *  *     ^ 


m,     itk+I  m  +1 

-1)    i81  93b1   L     C2+02b2(t1+l)] 


2     -i<t2<!       [2+e1b1(tl+i)+e2b2(t2+i)]  l 


rn,+2 


1 


max 


nu-Hiy-l    m.+2  ?    m,+2  m9+l    m9+l 

"1}  V     93bl        (V1)92        b2        C2+92b2(t2+l)] 


?ml+fll2  -l<t,<l 
-KtJ<l 


C2+81b1(t1+l)+82b2(t2+l)] 


m ,  +m  ?+4 


If    the    parameter    space,     9,     is     such     that       v,<8  .<*  v   <Q   <vi        and 


v3<93<w3>     then 
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max 
•Kt,<l 


|e(t1,t2) 


max 
2  -Kt,<l 


m„+l     m  +1 
w1w3b1(t1+l)w2        b2  £ 


2      -Kt^l  C2+v1b1(t1+l)+v2b2(t2+l): 


m?+2' 


wrJl+1vv3b^l+1[2+w2b2(t2+l)] 


2       Il<tJ<l  C2+v1b1(t1+l)+v2b2(t2+l)]   l 


m1+2l 


m,+2  2  m,+2 


nu+l  m?+l 


m,+m 


w!       w3bL       (t1+l)w2"     b26     [2+w2b2(t2+l)] 


1     2  -l<t,<l 


2   ^     "  -Kt^l  [2+v1b1(t1+l)+v2b2(t2+l)]   * 


m1-Hn?+4 


APPENDIX  2 
PROOF  OF  INEQUALITY  (2.31) 


Consider  model  (2.80)  for  which 

n(t;  _9)  =  a  +  (.49  -  a)exp(-S(15t  +  17))   ,       (A2.1) 

the  j-th  derivative  of     n(t;  _9)     with   respect  to  t   is 

n(j)(t;   9_)   =   (.49   -  a) ( -153)jexp (-B( 15t  +   17))      ,     j   =  2,3,4       . 

(A2.2) 

From   (A2.2)  we  have 


max    |n^J'(t;   8_)  |    =    |  .49   -  a|    |156|J   exp(-17S)     max    (exp(-15gt))    , 
•Kt<l  ~  -l<t<l 


j    =   2,3,4        .  (A2.3) 

Since  we  have   -l<t<l,   then  -K-t<l.   If  the  parameter  space,  9,  is 
such  that   .36<a<.41  and  .06<8<.16,  then  -150t<15S  and 
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93 


max   (exp(-15gt))   <;  exp(153] 
-Kt<l 


'A2.4) 


From  (A2.3)  and  (A2.4)  we  have 


max  |n 
•Kt<l 


;t;  6)|  <  (.49-a)(15)jeJexp(-26: 


2,3,4 


;A2.5) 


We  now  let  f(3)  =  8J'exp(-2S),  .06<S<.16,  j  =  2,3,4.  We  shall  find  an 
upper  bound  of  f(8)  over  the  interval  .06<3<.16.  The  first  derivative 
of   f(B)  with   respect   to   6   is 


f(1)(8)   =  jBJ'"1exp(-23)-2SJ'exp(-23)      ,     j    =  2,3,4 


We  set  f'^(S)  in  (A2.6)  equal  to  zero  and  solve  for  3  we  obtain 


A2.61 


8  -4 


j  -  2,3,4 


:A2.7] 


We  note  that  the  values  of  8  in  (A2.7)  are  equal  to  1,  1.5  and  2,  these 
values  maximize  f(3)  for  j  =  2,3,4,  respectively.  But  the  upper  bound 
of  8  is  .16.  Therefore,  8  =  .16  maximizes  f(3)  over  the  interval 
.06<S<.16  for  j  =  2,3,4.  We  substitute  S  =  .16  in  (A2.5)  and  obtain 


max  |nUj(t;  8)|<  (.49  -  a)  ( 15) j  (  .16)  jexp  (-  .32; 
•Kt<l 


j  =  2,3,4 


'A2.8) 
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Sy  substituting  (A2.8)  for  j  =  2,3,4  in  (2.66),  (2.67)  and  (2.63)  we 
obtain 


max  |g(t ;  9)  -  n(t ;  8) |  < 
-l<t<l 


;at)2(.49  -  a)(15)2(.16)2exp(-.32)/3  ,   for  linear  spline 


A2.9) 


<        (AT)3(.49-a)(15)3(.16)3exp(-.32)/8      ,       for  quadratic   spl i n< 


;at)4(.49   -  a)(15)4(.16)4exp(-.32)(5)/384     ,       for  cubic   spline 


Since    a    =    .36    maximizes     .49    -    a   over    the    interval       .36<a<.41,      then 
(A2.9)    becomes 


max      |g(t ;  _9)    -  n(t ;    8)  |    < 
•Kt<l 


;at)^(. 0679675)  for  linear  spline 


>t)  (.1631221197)  for  quadratic  spline 


(At)  (.0407805)  for  cubic  spline 


APPENDIX    3 
BOX  AND   LUCAS  DESIGN 


Consider  the  nonlinear  model  given  in  (2.1).  Box  and  Lucas's 
(1959)  design  criterion  for  estimating  _9  in  model  (2.1)  is  the  one  which 
maximizes  the  determinant  |MeL(_D,  _9q)  |  with  respect  to  the  design  matrix 
D_,  where  _e_Q  =  (91Q,  e2Q,  ...,  9  0)'  is  some  guessed  value  of  _9  and 
Mbi_(J1»  IflJ/o  is  Fisher's  information  matrix  of  order  pxp  with  Mn.  (D,  9) 
having  the  form 


n 
Mb.  (D,D    ■     I     h(x.,9)h'(x.,9)  (A3.1' 

i=1  1  -      -!    - 


where  x^  denotes  the  i -th  design  setting  of  x_(i=l,2,  ...,  n)  and  h_(x,9) 
is  the  partial  derivative  vector 


3f(x;  9] 
h(x,9)  = 


3e 


The    optimal     design     mentioned     above     is     said     to    be    D-optimal     in    the 
vicinity  of  9^   (locally   D-optimal). 
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